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FOREWORD 


This report was prepared by personnel of the Computational 
Mechanics Section of Lockheed’s Huntsville Engineering Center. 
It constitutes final documentation of efforts performed under 
Contract NAS8- 34507 for NASA-Marshall Space Flight Center. 

The NASA-MSFC Contracting Officer’s Representative for 
this research study was Dr. P.K. McConnaughey , ED32. 
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1. INTRODUCTION 

This report describes results of efforts by personnel of the Computational 
Mechanics Section at the Lockheed-Huntsville Engineering Center to assist the 
computational staff of NASA-MSFC in developing analytical capabilities for 
modeling hot gas flow on the fuel side of the SSME. The report specifically 
details work completed subsequent to the interim technical report submitted in 
August 1985 (Ref. 1). Computational and experimental results reported in Ref. 

1 will not be reiterated here, and the reader is referred to that earlier 
document for those details. 

Emphasis in this final report is placed on construction and documentation 
of a computational grid code for modeling an elliptical two-duct version of 
the Space Shuttle Main Engine (SSEM) fuel side hot gas manifold (HGM) . Also 
computational results for flow past a support strut in an annular channel. 

These three-dimensional results constitute the initial phase of a more detailed 
study of flow through the SSME/HGM strut region being completed by Lockheed 
under contract NAS8-37359 (Ref. 2). 

The approach for both of the aforementioned tasks is presented in Section 
2. Sample results are contained in Section 3. Included in Appendixes A, B, 
and C are a brief input guide for the two-duct HGM code, a listing of the code 
input file, and a source listing of the grid code itself. 
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2 . TECHNICAL APPROACH 


2 . 1 GENERAL 

The two-duct HGM grid generation uses as its fundamental building block 
the flexible algebraic techniques coded by the Lockheed Computational Mechanics 
Section into its own in-house geometry code. These techniques employ vector 
algebra, analytical geometry and transfinite interpolation to perform the 
computations in local curvilinear coordinates and generate the grid of discrete 
points in Cartesian coordinates. With this as a foundation, special sub- 
routines were constructed for describing particular surface shapes for this 
complicated structure such as the hold in the bowl and the fairing for the 
transfer duct. The combination of these codings was assembled, debugged, 
streamlined, and well commented for use by NASA-MSFC computational engineers. 

Early computational fluid dynamics results, reported in the interim 
report for this contract (Ref. 1) were performed with a finite difference code 
employing an explicit solution algorithm. Steady state was obtained from a 
trial initialization by performing successive iterations in time until all 
transients have involved away. The size of the time step used in this pro- 
cedure is a strong function of the density of points in the grid. The higher 
the density, the smaller the allowable time step. For this reason, relatively 
coarse grids were used, even for regions near the solid walls. 

Experience has dictated that much larger nodal densities near the walls 
are desirable for more accurate computational predictions. This precludes 
using an explicit code because of the unnecessarily large number of time steps 
required to obtain a steady state solution. The implicit code INS3D, developed 
at NASA-Ames (Ref. 3), was then chosen for additional SSME related computa- 
tions. The implicit solution algorithm incorporated into this code allows for 
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much larger time steps even for grids of large nodal densities. The results 
reported in Section 3.2 were obtained with INS3D. A brief description of this 
code, and how it was applied, is contained in Section 2.3. 


2.2 ELLIPTICAL TWO-DUCT MANIFOLD CODE 


Nearly all computational codes which are available to numerically solve 
the three-dimensional fluid flow equations are designed to be applied to a 
well structured grid model of the flow region. To perform the calculations, 
generalized independent variables are introduced Which transform the physical 
coordinates, (x, y, z), into general curvilinear coordinates, (n^, Hj, n^). 
Thus, the physical domain must be gridded as a single or series of hexahedral 
zones described by eight corner points, 12 edges, and six surfaces. Such an 
arbitrary zone is shown in Fig. 2-1. 



Fig. 2-1 Hexahedral Element Showing Local Intrinsic Coordinates 
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The approach used in the current study was to provide an algebraic grid 
generation code which would produce the Cartesian coordinates (x, y, z), for 
points along the lines of constant (n , n ). Basic mathematical 
techniques taken from analytic geometry and vector algebra were employed to 
describe a hexahedral zone in terms of piecewise continuous analytic functions 
which represent the zonal edges and surfaces. 

An intrinsic curvilinear coordinate system can be produced by mapping a 
unit cube onto the simply connected hexahedral zone. What is needed is a 
transformation function that will map a unit cube in (n^, space uni- 

valently onto the hexahedral volume of interest thus producing the required 
intrinsic coordinate system. A procedure which produces the desired result 
is referred to as either the method of transfinite interpolation of multi- 
variate blending function interpolation (Ref. 1). A brief description of 
this method follows. 


Let F (n, n., n,) be a vector-valued functional representing the region R 

X €m O 

of interest in curvilinear space. Then as <n , ti^ , n 3 > range over R, F traces 
out the region in Euclidean space (x, y, z) . Also let ♦ , ♦, and X be 
blending functions which obey the cardinality conditions: 


^(Tlj-1) 


♦j (n 2 =m) 

V’ , 3* n) 


j 1 if i = 1 
(o if i t 1 

j 1 if j = m 
( 0 if j ^ m 
( 1 if k » n 
(o if k ? n 


The simplest form of blending functions meeting these conditions are 


4> (t», ) * l - n. 

0 1 1 

vv 

vv 


l - n. 

I 

1 - n. 


♦ 1 (n ] ) = n i 

W = n 2 

W ■ " 3 


Then a trilinearly blended interpolant of F, which will map a unit cube onto 
the region R is given by 
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u(n lt n 2 ,n 3 ) = 


x(n 1 ,n 2 ,n ;j y 

y(n 1 ,n 2 ,n 3 ) 

z( WV' 


(l-n 1 )F(0 n 2 «n 3 ) + n 1 F(l,n 2 ,n 3 ) 


+ (l-n 2 )F(n 1 ,o,n 3 ) + n 2 F(n 1 ,i,n 3 ) 


+ (l-n 3 )F(n 1 ,Ti 2 ,o) + n 3 F(n 1> n 2 ,l) 


(i-n 1 )(l-n 2 ))F(0,o,n 3 ) - (i-n 1 )n 2 F(0,i,n 3 ) 


n 1 (i-n 2 )F(i,o,n 3 ) - n^Fd.i.^) 


(1 " n i )n 3 F<0,n 2 ,0) “ (l-n 1 )n 3 F(0,n 2 ,i) 


n 1 (i-ri 3 )F(l,n 2 ,0) - n 1 n 3 F(i,n 2 ,i) 


(i-n 2 )(i-n 3 )F(n 1 ,o,0) - (l-iyr^FO^o.i) 


- n 2 (i-n 3 )F(n 1 ,l,0) - n 2 n 3 F(n 1 .i.i) 


+ (i-n 1 )(i-n 2 )(i-n 3 )F(o,o,0) + (l-i^Xi-n^Fto.o.i) 


+ (i-n 1 )n 2 (i-n 3 )F(o,i.o) + (1 ~V V3 F(0,1,1) 


+ n 1 (i-n 2 )(i-ri 3 )F(i,o,0) + n 1 (l-n 2 )n 3 F(i,0,i) 


+ n n 2 <l-n 3 )F(l,i,o> + n 1 n 2 n 3 F(l,l,l) 


where 

F(0,n 2> n 3 ) represents a surface with 1)^=0, etc. 
F(0,0,n 3 ) represents an edge with n 1 *n 2 =0, etc. 
F(0,0,Q) represents a point with n 1 =n 2 *n 3 =0, etc. 
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In two dimensions this equation performs a bilinear interpolation over an 
arbitrary region consisting of four distinct corners simply connected by four 
edges, where 


F(0,n 2 ) = ty] 

along EDGE 4 

F(0 ,0) 


at P0INT 1 

F(i,n 2 ) = [y] 

along EDGE 2 

F(0,1) 

’ [ y' 

at POINT. 

A 

FO^.O) = ty] 

along EDGE^ 

F(1 ,0) 

■V 

at P0INT 2 

F(n 1 ,D = [yj 

along EDGE 3 

F(l,l) 


at P0INT 3 


and the interpolation equation for F could be rewritten as 

U = [*J « (l-n 1 )EDGE 4 + n 1 EDGE 2 + (l-n 2 )EDGE 1 + n 2 EDGE 3 

- (i-n 1 )(i-n 2 )POiNT 1 - d-n 1 )n 2 P0iMT 4 

- n 1 d-n 2 )P0iNT 2 - h 1 h 2 point 3 

In this equation EDGE, represents a vector-valued functional along edge four, 

4 

etc. Examination of this equation shows that it performs linear interpolations 

between EDGE, and EDGE„ and between EDGE., and EDGE . , hence the term bilinear 
13 2 4 

interpolation. Hence, if the functional for each edge can be derived and is 

analytic, a grid or mesh of node points can be generated by substituting values 

of n, and ti into the equation. 

12 

In three dimensions the general equation above performs trilinearly 
blended interpolation over an arbitrary region consisting of eight distinct 
corner points simply connected by 12 edges, where 
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F(0,n 2 ,n 3 ) * [^1 on SIDE 5 , etc. 

F(O f O,n 3 ) = [j[] along EDGE^ , etc. 

F(0 ,0,0) * [jf] at POINT^, etc. 
and which can be rewritten as 

U = tf] 

= (1-T) 1 )SIDE 5 + n 1 SIDE 6 + (l-n 2 )SIDE 2 + H 2 SIDE 4 + (1 -i> 3 )SIDE 1 + n 3 SIDE 3 

- (l-n i )(l-n 2 )EDGE 5 - (l-n i )n 2 EDGE g - U 1 (l-n 2 )EDGE 6 - n 1 n 2 EDGE ? 

- (i-n 1 )(i-n 3 )EDGE 4 - (i-n 1 )Ti 3 EDGE 12 - n 1 (i-n 3 )EDGE 2 - n 1 n 3 EDGE 10 

- (i-n 1 )(i-n 3 >EDGE 1 - (i-n 2 )n 3 EDGE 9 - n 2 (i-n 3 )EDGE 3 - n 2 n 3 EDGE u 

+ (1-n. ><i-n,) (l-n ) point + (l-n > (i-n,)n, point. + (l-n, )ti (l-n , > point 

L £ 3 1 1 4 j D 12 3 4 

♦ (i-n 1 )n 2 n 3 POiNT g + V 1_ V (1 *V P0Iirr 2 + n^i-n^n^oiin^) 

+ n 1 T) 2 (i-n 3 )POiNT 3 + n 1 n 2 n 3 POi*rr 7 
where SIDE., EDGE., POINT, represent vector-valued functionals 

1 J K 

on the surfaces, along the edges, and at the corner points, 
respectively. 
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This equation reduces to the previous two-dimensional analog along any 
flat surface or along any surface in which a straight line can be drawn 
between any two opposing edges such that the line lies entirely within the 
surface. 

With the general transformation, any point in local coordinates n , n^, 
can be related to the physical Cartesian coordinates x, y, z. The entire grid 
of discrete points is generated in the HGM code using this concept. This 
general interpolant can accommodate any stretching function for concentrating 
points near walls or regions of large gradients. Furthermore, the edges of the 
hexahedral can be segmented allowing another means of grid spacing control. 

2.3 STRUT IN ANNULUS 

Previous SSME/HGM computations reported earlier by Lockheed (Refs. 1 and 
4) and other investigations (Ref. 5) have either not modeled the support strut 
region in the manifold or poorly approximated its influence on the flowfield 
environment. It was for this reason that NASA requested an additional compu- 
tations be made which accurately included these obstacles in the flow path at 
the entrance to the fuel bowl of the manifold. Initial work on this task was 
performed under this contract and is being reported here. Results of follow-on 
work are presented in the final report for NAS8-37359 (Ref. 2). 

The approach was to generate a computational grid consisting of a single 
strut in an annular channel with the dimensions of the strut size, channel 
width, and channel curvature being approximately the same as those in the 
actual SSME/HGM. Numerical experiments were then to be performed using a 
three-dimensional incompressible Navier-Stokes code which employed an implicit 
solution algorithm. The Lockheed in-house algebraic grid code was used to 
model the geometry and INS3D was used to obtain flowfield solutions. 
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The INS3D code solves the three-dimensional incompressible Navier-Stokes 
equations in primitive variables. An implicit finite difference operator is 
used in a general curvilinear coordinate system. The solution procedure uses 
the standard approximate factorization scheme. The pressure field solution is 
based on the concept of adding a time-like pressure term into the continuity 
equation via an artificial compressibility factor. This approach was first 
introduced by Chorin (Ref. 6) and later adopted by Steger and Kutler (Ref. 7) 
using an implicit approximate factorization scheme by Beam and Warming (Ref. 

8). It is from these earlier developments that INS3D evolved (Ref. 3). 

Values of the artificial compressibility factor are bounded in order not 
to influence the steady state mass conservation. In the INS3D methodology 
mass conservation is of crucial importance if a stable solution is to result. 
Since the continuity equation is modified to obtain a hyperbolic- type equation, 
pressure waves of finite speed will be introduced. The speed of propagation 
of these pressure waves depends on the magnitude of the compressibility 
parameter. When the pressure waves travel through a given location a pressure 
gradient is created there. Near boundaries, the viscous boundary layer must 
respond to this pressure fluctuation. To accelerate convergence and avoid 
slow fluctuations it is desirable that the time required for pressure waves to 
propagate through the region of interest be much less than the time needed for 
the boundary layer to fully adjust itself. This condition provides for a 
lower bound on the artificial compressibility factor. The upper bound on this 
factor comes not from the physics but from the effects of the approximate 
factorization of the governing equations. When the finite difference form of 
the equation is factored, higher order cross-differencing terms are added to 
the left-hand side of the equation. These added terms must be made smaller 
than the original terms everywhere in the computational domain. This condition 
results in an upper bound on the compressibility factor. 
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It is well know in the computational fluid dynamics community that the 
approximate factorization schemes which employ alternating direction type 
implicit methods have stability problems in three dimensions (Refs. 9 through 
12). The INS3D code satisfactorily overcomes this difficulty by providing 
second and fourth order smoothing terms to the algorithm to ensure stability 
without adversely affecting mass conservation. 

Currently, the Computational Mechanics Section at Lockheed-Huntsville has 
the INS3D code operational on the Cray-XMP at NASA-MSFC and NASA-Ames as well 
as on its own VAX 11/785. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


3. RESULTS 


3.1 MANIFOLD GRID CODE 

3.1.1 Geometry and Grid 

A schematic representation showing the construction of the geometry which 
has been modeled is presented in Fig. 3-1. Only the outer wall is displayed and 
two perspectives are shown. These solid surface plots were generated with the 
actual output from the grid code. Only half of the manifold is represented 
since, due to the plane of symmetry Which divides the two transfer ducts, only 
half need be computationally modeled. 



Fig. 3-1 Schematic Showing Outer Wall of HGM Geometry that has been Modeled 
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The geometry is generated in five separate pieces or zones. The five 
zones are shown in Figs. 3-2 and 3-3. Each zone has a general hexahedral 
shape, and in these figures all of the eight corner points of each are clearly 
indicated. Figure 3-4 shows the Cartesian coordinate system relative to which 
the position of each node in the grid is referenced. Except for the hole 
perimeters in the outer wall of Zone 2 and the fairing at the entrances to 
Zone 3, the edges of each zone can be described as piecewise continuous seg- 
ments composed of either straight lines, circular arcs, or elliptical arcs. 

In addition, excluding the two special regions previously mentioned, the sur- 
faces of each zone can be generated by rotating an edge about an appropriate 
axis. For Zones 1, 2, 4, and 5, the X axis is the axis of revolution. 

Zones 1 and 2 make up the bowl section of the manifold and contain 58 
nodes in the x-direction, from bowl entrance to rear of the bowl, 109 nodes in 
the circumferential direction and 21 nodes between the inner and outer sur- 
faces. Distribution of nodes in the bowl is presented in Fig. 3-5. 

Zones 4 and 5 , the turnaround duct (TAD) , are composed of 71 nodes in the 
streamwise direction, 109 nodes in the circumferential direction (0 to 180 
deg), and 21 nodes across the duct between inner and outer surfaces. Two 
perspectives showing this nodal distribution are shown in Fig. 3-6. 

Zone 3 is the elliptical transfer duct portion. This duct has been gen- 
erated with 59 nodes along the duct axis and 44 x 30 nodes in a cross section. 
Surface and cross-section grids for this zone are given in Fig. 3-7. Figure 
3-8 displays a view of the manifold outer wall grid to show grid continuity 
from one zone to another. To summarize, each zone contains the following 
number of nodal points: 

e Zones 1 and 2 (Bowl): 132,762 

e Zones 3 (Transfer Duct): 77,880 

e Zones 4 and 5 (TAD): 162,519 

e Manifold total Nodes: 373,161 
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Fig. 3-2 Schematic Showing First Two Zones (Zone 1 Left; Zone 2, Right) 
into Which the Two-Duct HGM was Subdivided 
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Fig. 3-3 Schematic Showing Zones 3, Top, and A, Middle, 
and 5 of the Two-Duct HGM Model 
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Fig. 3-4 Cartesian Coordinate System Relative to 

Which Positions of All Modes are Referred 
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Fig. 3-5 Grid Plots for Inner and Outer Surfaces of Zones 1 and 2 

as Well as Internal Grid at the Common Place of Intersection 
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Fig. 3-6 Grid Plots of TAD Internal Grid (Bottom) and Grid 
Distribution on Inner and Outer Surfaces 
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Fig. 3-7 Surface and Cross-Section Node Distribution 
for Transfer Duct 
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Fig. 3-8 Grid Plot Showing Node Distribution on Outer Wall of Zones 1, 2, 

and 3 Displaying Continuity of Grid Lines from One Zone to Another 
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The internal grid resolution presented in the previous figures is adequate 
for a laminar viscous computation but would need to be modified for application 
to a turbulent computation. The procedure for doing this will be described in 
Section 3.1.4. 

3.1.2 Computer Code 

A concise input guide and an input listing for the grids displayed in 
Figs. 3-4 through 3-8 is also included in Appendixes A and B, respectively. 

A source listing of the two-duct HGM geometry code is provided in Appendix C. 


Figure 3-9 shows the primary calling sequence. For completeness, 
secondary calling sequences are given in Fig. 3-10. A brief explanation of 
the function of the major subroutines in the primary calling sequence follows: 


e INITIAL Reads the first two lines of the input file, initializes 
coefficient arrays, and defines logical unit numbers and 
counters . 

e INPUT Reads the remainder of the input file and sets all parameters 
to be used in remaining subroutines. 

e MESH The controlling subroutine for the generation of the spatial 
coordinates of each node in each zone. 


e ETABC Calculates the values of nj., 1 ) 2 * **3 al° n S the I, J, and K 
directions for each hexahedral shaped section of each zone. 

e EDGES Determines the Cartesian coordinates for nodes along each edge 
of each side of each section of each zone using the bilinear/ 
trilinear interpolation scheme. 

e SURFACE Determines the Cartesian coordinates of nodes on a three- 

dimensional surface using the trilinear interpolation scheme. 
Here, all outer and internal surface nodes are calculated from 
the previously determined edge distributions. 


e OUTPUT Provides printed output and stores geometry in File 20 for use 
in plotting or as input to an integration code. 


The output to File 20 is in the format to be input as a multi-grid 
geometry file to the PL0T3D plotting code. 
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Fig. 3-10 Secondary Calling' Sequences for Two-Duct HGM Grid Code 


3-12 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 






LMSC-HEC TR F225968 


The grid code is written in standard FORTRAN. However, the listing pro- 
vided in Appendix C is a VAX 11/785 version and could contain some generic VAX 
statements which would have to be translated if used on a different machine. 

3.1.3 Code Implementation 

To implement the HGM geometry code as it appears in the listing of 
Appendix C, the input listing of Appendix B must correspond to logical Unit 5. 
Unit 6 must be assigned to the written output and Unit 20 to the geometry file 
Which will contain the x, y, z coordinates for each node in each zone. 

The input file is labeled on the card image which begins each zone input 
to indicate the zone being described by the succeeding card images. These 
zone labels correspond to those shown in Figs. 3-2 and 3-3. The input file, 
as listed in Appendix B, will create a multi-grid file containing all five 
zones. Using this input is, of course, possible only if run on a computer 
with sufficient CPU and/or disk storage. Alternatively the bowl, transfer 
duct, and TAD can be run separately by changing the second parameter on card 
image two, of the input, from a 6 to a 1, 3, and 4, respectively. 

A detailed description of the input file to the code is provided in 
Appendix A. Modifications to the geometry can be facilitated by studying the 
input guide while observing both the input file and the detailed grid pictures 
presented in Figs. 3-3 through 3-8. Redistributions of nodes can be accomp- 
lished by making minor modifications to the input files. For example, for a 
turbulent computation if the nodes in Zone 4 near the wall require redistribu- 
tion closer to the wall then a change would need to be made to card type 9 on 
the last line of the input for that zone. The 6.0 appearing in the n 2 
position could be changed to 10.0 (see page A-10). 

Dimensioning in the program has been kept to a minimum. The largest 
dimensioned array in the bulk of the code is NODENUM(5000) . At the very end 
of all computations the PL0T3D file can be generated by employing the program 
in Appendix D. 
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In this subroutine the x, y, z coordinate arrays are each dimensioned to 

200,000. The 5,000 corresponds to twice the maximum number of nodes in a 

plane perpendicular to the marching direction (q direction input on card 7) 

for creating the geometry. The 200,000 must be equal to or greater than the 

number of nodes in the largest zone, which is zone 5, the second half of the 

TAD. In Zones 1, 2, 4, and 5, the direction of is from TAD entrance to bowl 

back wall; the ti direction is from inner to outer wall; and ti is directed 
2 3 

from side opposite transfer duct circumferentially. In zone 3, increases 
in the streamwise direction from bowl outward, and x form the cross 
planes in the duct. 

Note that the code is designed to output each zone of the manifold so that 
each has one cross plane in common with the preceding zone(s) . This must be 
remembered for incorporating the grid into a flowfield solver code. The 
geometry must be integrated in a multi-block or multi-zone fashion. If the 
computer available has large enough core memory or if it is a large virtual 
machine then Zones 1 and 2, and 4 and 5 can easily be combined into two larger 
zones since at all common planes the nodal positions match exactly. It is not 
recommended that zones 1, 2, 4, and 5 be combined into one zone. This could 
easily be accomplished, since continuity in all three n directions would be 
maintained, but the total number of nodes would be untenable on all but a Cray 
II or an ETA 10 machine. 

3.2 STRUT LAMINAR COMPUTATION 


The Lockheed-Huntsville algebraic computational grid generation code was 
employed to generate a model for a single support strut in an annular channel. 
A C-type grid was selected for nodal distribution in parallel annular surfaces 
and is shown in Fig. 3-11. This is a 26 x 201 node structure with stretching 
toward the strut surface. The three-dimensional model consisted of 31 such 
surfaces in concentric cylindrical fashion with stretching toward inner and 
outer surfaces. Figure 3-11 also shows part of the solid surfaces in the 
model. The total number of grid points in the computational domain was 
162,006. 
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Fig. 3-11 Grid Plots Showing C-Grid Used in Each of the 31 
Circumferential Planes (Top) and Partial Surface 
Grids of Strut and Inner Outer Wall 


3-15 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 




LMSC-HEC TR F225968 


The basic INS3D code was modified for treating the sides of the compu- 
tational domain as periodic boundaries and for special treatment of the common 
grid plane at the center of the grid behind the strut. At the entrance plane, 
the velocity components were initialized for fully developed laminar channel 
flow at zero degree incidence to the strut. Velocity and pressure were then 
held fixed throughout the computation. A nondimensional time step of 0.025 
was used and a steady solution for Reynolds number of 500 was obtained in 800 
iterations. Figures 3-12 to 3-21 present various characteristics of this 
steady laminar solution. Velocity magnitude contours and static pressure con- 
tours in the region surrounding the strut in concentric cylindrical surfaces 
near the inner wall, center of annulus, and near outer wall are shown in Figs. 
3-12, 3-13, and 3-14, respectively. An expanded view of the same information 
for the central annular surface appears in Fig. 3-15. The velocity contours 
in the figure clearly shows the extension of the strut wake region several 
strut lengths downstream. In the SSME/HGM the support struts are at the 
entrance to the bowl (exit of the TAD) . Even though these current results are 
laminar and steady, the influence of the strut wake is shown to be significant 
and could cause a significant difference in the predicted flow through the 
transfer ducts. 

More details of the strut near wake region for a central circumferential 
surface and central radial surface are given in Figs. 3-16 and 3.17. The in- 
fluence of the wake flow pattern on flow particles originating upstream of the 
wake is indicated by the particle traces shown in Fig. 3-18. To further illus 
trate the three-dimensional character of this region of the wake Fig. 3-19 and 
3-20 trace particle paths beginning at positions in the wake itself. Clearly, 
if the wake were unsteady, as it is in the actual HGM, these complicated flow 
patterns will move into the bowl, interact, and exit through the transfer 
ducts. 

For completeness, surface pressure distribution on the strut itself is 
presented in Fig. 3-21. 
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Fig. 3-12 Static Pressure Contours (Top) and Velocity Magnitude 
Contours for Strut in Annular Channel near Inner Wall 
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Fig. 3-13 Static Pressure Contours (Top) and Velocity Magnitude 
Contours for Strut in Annular Channel at Mid-Channel 
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Fig. 3-14 Static Pressure Contours (Top) and Velocity Magnitude 
Contours for Strut in Annular Channel near Upper Wall 
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Fig. 3-15 Static Pressure Contours (Top) and Velocity Magnitude 

Contours for Entire Length of Computational Circumferential 
Plane at Mid-Channel 
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Velocity Vectors in near Wake Region in a Central Plane for Case 
of Strut' in Annular Channel (Mean Flow is from Left to Right) 




Fig. 3-17 Velocity Vectors (Top) and Velocity Magnitude Contours 

in Central Radial Plane in near Wake Region Behind Strut 
in Annular Channel (Mean Flow is from Left to Right) 
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Fig. 3-18 Particle Traces for Flow Past Strut in Annular Channel 
Showing Paths of Particles Released near Front of Strut 
in a Plane near the Inner Wall (Top) , one Quarter Channel 
Height Above Inner Wall (Bottom Left) and Mid-Channel 
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Fig. 3-19 Particle Traces of Particles Released near Rear Surface of Strut 
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Fig. 3-20 Particle Traces of Particles Released near Central Radial 
Plane but just Downstream of Strut Rear Surface 
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Fig. 3-21 Surface Pressure Contours on Strut Front Surface 
(Left, Flow into Page) and Side Surface (Right, 
Flow Left to Right) 
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4 . CONCLUDING REMARKS 

An elliptical two-duct SSME fuel side hot gas manifold geometry code has 
been developed for use by the Computational Fluid Dynamics Staff of NASA-MSFC. 
This report describes the methodology of the program, makes recommendations on 
the implementation, and provides an input guide, input deck listing, and source 
code listing for the code. The listing is well commented in order to assist 
the user in following its development and logic. A magnetic tape containing 
the source deck will be provided, upon request, to NASA-MSFC for use on its 
EADS network. 

The NASA-Ames three-dimensional incompressible Navier-Stokes computa- 
tional fluid dynamics code, INS3D, was obtained and implemented on the MSFC 
IBM/Cray computer facility. A low Reynolds number laminar calculation was 
performed with this code on a 162,006 node model of a strut in an annular 
channel. The dimensions were approximately those of the SSME/HGM fuel bowl 
entrance region which contains 12 such support struts circumferentially 
distributed in hot gas flow path. The computation was made as an initial step 
in a thorough numerical investigation of the influence of these obstacles on 
the flow exiting the manifold and impinging on the main injector LOX posts. 

Results of this steady laminar computation indicate that a complete 
three-dimensional analysis of the whole manifold would require a "strut zone" 
for reliable predictions of of the transfer duct exit plane flow structure. 
This is especially evident since flow visualization results have shown that 
the duct flow is largely unsteady (Ref. 13). 
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Appendix A 


INTRODUCTION 

This geometry input guide is presented in two sections: (1) a definition 
of terminology commonly used for inputting and describing the geometry, and (2) 
a summary of card types used to input the geometry and a detailed description 
of the associated parameters and their input values. 

We begin with an overview of how to apply the program. The flowfield 
domain is divided into zones in order to simplify the input necessary to 
describe the complicated geometry. Each zone contains its own internal 
coordinate system, and is described using points, edges, and surfaces. An 
edge may consist of from one to ten segments. A segment or a surface may 
require special input depending on its type. This is the case for zones 2 and 
3 of this two-duct HGM model. 

The second section of this appendix presents a detailed description of 
the input parameters. Each card type is listed in the order of input with its 
associated parameters. Each parameter is identified as to its usage in the 
program with the options of each shown. Reference to Fig. A-l or Table A-l 
may be necessary to explain some of the input parameters and their order of 
input. All of the card types are not necessarily input for a specific zone, 
but may be set in specific subroutines in the program itself. 

Card type 9 may be used when other than an equal distribution of nodes is 
desired, etc. Whereas card type 13 is necessary if there will be more than one 
segment per edge. And cards type 10 and 11 are used if additional information 
is needed to describe a segment or a surface. Certain of the input parameters 
on early cards dictate which of the later cards are read in. 
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DEFINITION OF TERMINOLOGY 


Edge 


An edge consists of from one to five segments. Twelve edges are 
used to describe a 3D zone. 


Map The geometry maps a point from n space into real space. When 

describing a surface mapping one could say that setting map = 2 
refers to a planar n space surface being mapped onto a cylin- 
drical surface in real space. 


Node At each intersection of n coordinates a node is generated by 

the program forming the grid which will describe the flowfield 
domain . 


Point 



The comers of a section are called ’‘points.” The location and 
initial flow directions are input for each point. There are 
eight points. 


Segment An edge is subdivided into as many as five segments. A segment 

may be a straight line, a circular arc, a helical coil, a 
trigonometric function of angle or length, a cubic spline, or 
user defined. 


Surface A three-dimensional section will consist of six surfaces which 

form a generalized hexahedron. A surface may be planar, 
cylindrical, an edge of revolution, or user defined. 


Zone The flowfield domain may be subdivided into zones. Zones are 

generated independently and are the fundamental building block. 
Each zone contains its own n coordinate system. 


n^, ^2 > H3 Localized coordinate directions within a zone. These co- 
ordinates describe a cube in T) space. 0 < t)i < 1 
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SUMMARY OF CARDS 


Card Type Parameter List/Format 


Flowfield 

Parameters 


1 ITITLE(I), 1=1,80 

( 20A4 ) 

2 NZONE, IZINDEX, MAPTEN, INCHES 

(815) 

3 DCTl(I), DCT2 ( I ) , DCT3(I), 1=1,37 

(8F10.3) 


Zone Parameter 


4 

Section 

Parameters 

NSECT 

(15) 


5 

MAPEDGE(I), 1=1,12 

(1215) or (3(4110)) 


6 

MAPSIDE(I) , 1=1,6 
(615) 


7 

MARCH 

(15) 


8 

(NMBRNDS(I), 1=1,3), (ISTRTCH(I), 1=1,3) 
(615) 


9 

STRETCH ( I ) , 1=1,3 
(3E10.4) 


10 

[ (COEFE(I ,K, J) , 1=1,8), K=1 , 5 ] , J=1,4(2D) OR 
(8E10.4) 

12(3D) 

11 

[COEFS(I, J) ,1=1,8] , J=l,6 
(8E10.4) 


12 

(POINT(I, J) , 1=1,3) ,J=1, 8 
(8E10.4) 


13 

[ (SEGMAXd ,K, J) , 1=1,3), ETAMAX(K.J), K=1,4J 
(6E10.4) 

J=l, 12 
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CARD AND VARIABLE DESCRIPTIONS 


Flowfield Parameters 


CARD TYPE 1 Problem Identification Label 


Format (20A4) 


ITITLE Alphanumeric information used for identifying the flowfield 
geometry. Columns 1-80 are read and printed only. 


CARD TYPE 2 Problem Option Controls Flags 


Format(5I5) 


NZONE The number of zones into Which the flowfield geometry is 
divided. The maximum number of zones is 99. 

izindex Zone index for selecting individual components of the manifold 
to be computed. 

= 1 zones 1 and 2 of Bowl 
= 3, Transfer Duct 

= 4 zones 4 and 5 of Turnaround Duct 
= 6 all five zones 

MAPTEN This option determines the maximum number of segments which 
will be input per edge. 

= 0 Five segments per edge Format (1215) 

= 1 Ten segments per edge Format(3(4I10) ) 

INCHES This option specifies the dimensions of the coordinates being 
input. The output data will be written in feet for compati- 
bility with the INTEGRATION program. 


a o Dimensions in feet 
a 1 Dimensions in inches 
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Bowl Hole Parameters 


CARDS TYPE 3a Angles, in degrees, about the axis of the transfer duct, 

measured counterclockwise from the x, z-plane as viewed down 
the z-axis. 


CARDS TYPE 3b Radial distance difference (in inches) between transfer duct 
ellipse and hole perimeter along angular directions specified 
on cards 3a. 

CARDS TYPE 3c Radius of curvature (in inches) at each angular distance, 

specified on cards 3a, for describing the transfer duct weld or 
fairing. 


Zone Parameter 


CARD TYPE 4 Index dividing zone data Format (15) 

NSECT Integer must be 1 and is used for separating data for each zone 
on the input file. 


Section Parameters 


CARD TYPE 5 Edge Shape Function Indicators Format(12I5) or 3(110) 

HAPEDGE(I) . 1=1. 12 

These are packed integer flags that specify which edge shape 
functions will be used for the current section. The edges are 
input in numerical order. The edge numbers are defined 
according to Fig. A-l. The user should study this figure 
before inputting the geometry. 

Each of the edges may consist of up to ten segments with each 
of these segments having its own shape function. The value of 
MAPEDGE(I) can consist of up to ten integers packed into one 
word MAPEDGE(I). HAPTEN specifies the maximum number of seg- 
ments per edge. The edge shape function indicators for each 
segment are input in chronological order of increasing tj for 
each edge with the final packed integer being right adjusted. 
For example, if HAPEDGE(4) = 112, then edge 4 consists of three 
segments: the first segment is type 1; the second segment is 
type 1; and the third segment is type 2. If only one segment 
describes an edge, then only one indicator is used, right 
adjusted. 
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A library of edge shape functions indicators for the HGM 
GEOMETRY program follow. If any edge shape function other than 
a linear segment is specified, then edge coefficients (COEFE(I)) 
must be input. Card type 1 is used to define the analytical 
function describing a segment. 


= 1 

Linear segment 



= 2 

Circular arc 

( input 

COEFE(I) ) 

= 3 

Conics 

( input 

COEFE(I)) 

= 4 

Edge of revolution 

( input 

COEFE(I) ) 

= 5 

Special segment 

( input 

COEFE(I)) 

= 6 

Special segment 

( input 

COEFE(I) ) 

= 7 

Special segment 

( input 

COEFE(I)) 


CARD TYPE 6 Surface Shape Function Indicators Format(6I5) 

MAPSIDE(I) .1=1.6 

These are integer flags that specify which surface shape 
functions will be used for the current section. These flags 
are input only for three-dimensional problems since two- 
dimensional geometries are defined completely by the edge 
functions. The surfaces are input in numerical order. The 
surface numbers are defined in Fig. A-l. The user should study 
this figure before inputting the geometry. An edge of revolu- 
tion requires the input of surface coefficients (COEFS(I)) on 
card type 15 to define a relative origin on the axis, the axis 
of revolution, and the direction of revolution. 

= 1 Planar surface 
= 2 Cylindrical surface 

= 3 Special surface (user defined) 

= 4 Edge of revolution (input COEFS(I)) 

=5 Hole in bowl surface 
= 6 Duct surface at bowl 
= 7 Duct fairing surface 
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CARD TYPE 7 Node Numbering Sequence Specs Format(6I5) 

MARCH (15) The value of MARCH determines the node generation and hence the 
node numbering sequence, (default = 1) 

The numbering sequence corresponding to follows. 



= i 

TI3, Ti2» n l (default) 



= 2 

hit h3» t)2 



= 3 

t > 2» ”3 


CARD TYPE 8 

Node 

Distribution Parameters 

Format(6I5) 


NMBRNDS ( I ) . 1=1 3 

NMBRNDS ( I ) is the number of nodes in the m direction for 
the current section. The limit is 200 nodes in any coordinate 
direction. This may be changed in the program by respecifying 
the ETAS (3, 200) array. 

ISTRTCH(I) .1=1 3 

This option gives the user control over the node distribution 
in each of the coordinate directions. 

= 0 Uniform spacing 

= 1 Input actual nx values for NMBRNDS(I) nodes 
(input ETAS ( I ) ) 

= 2 Decrease spacing in iij direction. Input a stretching 
factor greater than 0.0 in STRETCH(I) . 

= 3 Increase spacing in t)j direction. Input a stretching 
factor greater than 0.0 in STRETCH(I). 

= 4 Double stretching. Input a stretching factor greater than 
0.0 in STRETCH(I). Use an odd number of nodes. 

= 5 Decrease spacing in t»i direction. Input minimum grid 

spacing as a percentage of the total length in STRETCH(I) . 
= 6 Increase spacing in tix direction. Input minimum grid 

spacing as a percentage of the total length in STRETCH(I). 
= 7 Double stretching. Input minimum grid spacing as a 

percentage of the total length in STRETCH(I). Use an odd 
number of nodes. 

If ISTRTCH(I) * 1, input a set of cards type 13 for each q 

direction to be input. 

If ISTRTCH(I) > 2, input card type 12. 
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CARD TYPE 9 Option for Stretching Function Format(3E10. A) 

(input when ISTRTCH(I) > 2) 

STRETCH (I). 1=1. 3 

This parameter is input for each coordinate direction 
designated for stretching by ISTRTCH(l) > 2. 

Example: Several stretching functions will be demonstrated using 21 

points for comparison. Note, that total length * 10.0 for 
ISTRTCH = 6 and 7. 

ISTRTCH = 3 STRETCH » 2.0 


STRETCH =4.0 


STRETCH * 6.0 



I 
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STRETCH * 4.0 


STRETCH « 


STRETCH > 8.0 


STRETCH - 10.0 



ISTRTCH » 6 STRETCH • .002 


m il l m M i., 

STRETCH « .004 


STRETCH - .004 


STRETCH a .008 
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STRETCH - .010 


ISTRTCH « 7 


STRETCH • .002 


STRETCH ■ .004 


STRETCH - .006 


STRETCH - .008 


STRETCH - .010 


CARD TYPE 14 Coefficients for Edge Shape Functions Format (8E10. 4) 

(input for each segment in MAPEDGE(I) > 1) 

COEFE(I) .I=»1.8 

These coefficients are used to describe the edge shape 
functions for each segment of the current section. The 
coefficients for each segment are input on separate cards in 
the same order as the indicators on card type 5. 
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Map Type Input Parameters 

2 Circular Arc COEFE(I), 1=1,3 are the x, y, and z coordinates of the 
center of the arc. 


3 Edge of COEFE(I) , 1=1,3 are coordinates of the center 

Revolution of the arc. 

COEFE(I) , 1*4,6 are components of the unit vector 
along the axis. Direction according to right hand 
rule when a vector is revolved from point 1 to point 2. 


CARD TYPE 15 Coefficients for Surface Shape Functions Format (8E10. 4) 
(input when MAPSIDE(I) = 4 and IDIM * 3) 

COEFS(I) .1=1.8 

These are the coefficients defining the surface shape functions 
for each surface formed by an edge of revolution in the current 
section. Each surface which has MAPSIDE(I) = 4 on card type 8 
is input on a separate card in the same order as they occur on 


Input Parameters 

Surface formed by revolving an edge about an 
axis. 

COEFS(I), 1*1,3 a point on the axis of revolu- 
tion which becomes the origin of a local co- 
ordinate system. This point must lie outside 
of the projection of the edge onto the axis of 
revolution. 

COEFS(I) , 1*4,6 are components of the unit 
vector along the axis of revolution in the 
direction of increasing n. 

COEFS ( 7 ) indicates the rij direction in 
Which the edge will revolve. 


CARD TYPE 16 Coordinates of Points Format (5E10. 4) 

POINT(I.J). J * Point number 

These parameters are the coordinates and flow direction at each 
comer of a general hexahedral(3D) . Figure A-l shows this con- 
figuration with the points numbered from 1 to 8. There are 
eight cards of type 16 to be input. 


card type 8. 
Map Type 
4 


Surface of 
Revolution 
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POINT ( 1 , J ) - the x coordinate of point J 

POINT(2 , J) - the y coordinate of point J 

POINT ( 3 , J ) - the z coordinate of point J 

Important Note: All cards type 12 are not input consecutively. 

They are grouped with cards type 13. See Table A-l for the 
exact sequence of card types 11 and 13. 


CARD TYPE 17 Segment Extremals for Edges Format(6E10.4) 

SEGMAXd.K. J) . K = Segment Number. J = Edge Number 

Each edge may be segmented up to five times. Therefore, cards 
type 13 are repeated for each successive segment of edge J. 

Each segment must be input on a separate card type 13. The 
extremal for the final segment of an edge is not to be input 
since this point is already defined by the POINT(I) input. The 
number of cards type 13 for each edge will thus be one less 
than the number of segments on that edge. In particular, if an 
edge consists of only one segment, no cards of type 13 are 
input for that edge. 

See Table 1 for the input order of card types 12 and 13. Each 
POINT(I) is input on a single card, followed by up to five 
cards containing the extremals. 

SEGMAX(1,K, J) - The extremal x coordinate for the 
segment of edge J. 

SEGMAX(2 ,K, J) - The extremal y coordinate for the K^ 1 
segment of edge J. 

SEGMAX ( 3 , K , J ) - The extremal z coordinate for the K th 
segment of edge J. 

ETAMAKK. J) . K= Segment Number. J= Edge Number 

The maximum value of the nj coordinate on the K tli segment 
of edge J. Input a negative value when defining a fold line. 
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TWO DUCT HGM 11+ 


2 

i 

0 

i 









0 

.0 

10. 

0 

20 

.0 

30.0 

40.0 

50.0 

60.0 


70.0 

80 

.0 

90. 

0 

100 

.0 

110.0 

120.0 

130.0 

140.0 


150.0 

160 

.0 

170. 

0 

180 

.0 

190.0 

200.0 

210.0 

220.0 


230.0 

240 

.0 

250. 

0 

260 

.0 

270.0 

280.0 

290.0 

298.0 


310.0 

320 

.0 

330. 

0 

340 

.0 

350.0 

360.0 

ANGLE 




0 

.30 

0. 

29 

0 

.28 

0.25 

0.23 

0.20 

0.16 


0.12 

0 

.06 

0. 

05 

0 

.06 

0.12 

0.20 

0.28 

0.37 


0.46 

0 

.55 

0. 

66 

0 

.76 

0.90 

1.25 

1.85 

2.45 


2.35 

1 

.85 

1. 

32 

1 

.00 

0.95 

1.05 

1.30 

1.649 


1.00 

0 

.75 

0. 

55 

0 

.43 

0.35 

0.30 

DIFFERENCE 



2 

.0 

2. 

0 

2 

.0 

2.0 

2.0 

2.0 

2.0 


2.0 

2 

.0 

2. 

0 

2 

.0 

2.535 

3.069 

3.604 

4.139 


4.673 

4 

.673 

4. 

673 

4 

.673 

3.81 

3.151 

2.744 

2.684 


2.801 

3 

.055 

3. 

596 

4 

.2 

4.253 

3.509 

2.54 

1.942 


2.0 

2 

.0 

2. 

0 

2 

.0 

2.0 

2.0 

CURVATURE 



1 








BOWL 

ZONE 

i 


111 

1 

113 

1 

4 

4 

4 4 

111 111333 1 




1 

•j 

4 

1 

4 

1 

4 







58 

21 

36 

0 

4 

0 








0.0 


6.0 


0.0 







1 







BOWL WITH 

HOLE 

ZONE 

2 


111 

111333 

1 

4 

4 

444 444 

111 111113 1 




1 

4 

1 

5 

1 

1 







j 

58 

21 

74 

0 

4 

0 








0.0 


6.0 


0.0 







1 








DUCT 

ZONE 

3 


7 

91 

8 

91 

7 

7 

8 8 

7 91 

8 91 




7 

9 

6 

7 

1 

7 

7 







30 

59 

44 

4 

0 

4 








6.0 


0.0 


4.0 







1 








TAD t 1 

ZONE 

4 


12 

1 

12 

1 

22 

22 

2222 2222 

12 1 

12 1 




1 

2 

1 

2 

1 

2 







1 












30 

21 

109 

0 

0 

0 







1 








TAD t 2 

ZONE 

5 


211 

1 

211 

1 

22 

22 

2222 2222 

211 1 

211 1 




1 

2 

1 

2 

2 

1 







i 

42 

21 

109 

0 

0 

0 








/EOR 
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****MAIN 


PROGRAM HGM2DUCT( INPUT, OUTPUT, FILE20, TAPES, TAP E 6 -OUTPUT , 
& TAPE20-FILE20) 

PROGRAM HGM2DUCT 


ELLIPTICAL TWO-DUCT HOT GAS MANIFOLD GRID CODE 
DEVELOPED BY THE COMPUTATIONAL MECHANICS SECTION 
LOCKHEED ENGINEERING CENTER, HUNTSVILLE, ALABAMA. 
TAPE20 GEOMETRY DATA 


COMMON /COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
COMMON /INITA/ I Z INDEX, MAPTEN, INCHES 
COMMON /ZONING/ I ZONE , ISECT,NZ INDEX, NMBRNDS< 3 ) 
COMMON /UNITS/ NU5,NU6,NU20 

— INITIALIZE PROGRAM 

CALL INITIAL(NZONE) 


ZONE 

DO 300 IZONE-l,NZONE 
READ( NU5 , 1000 ) NZINDEX 

READ INPUT FOR EACH ZONE 

CALL INPUT 

GENERATE MESH FOR EACH ZONE 

CALL MESH 
300 CONTINUE 

PRINT FILE STATUS 

CALL STATUS ( 0 ) 

STOP 

FORMAT STATEMENTS 

1000 FORMAT (15) 

END 

A********************************************************** 

***** INPUT**** ********************************************* 

SUBROUTINE INITIAL(NZONE) 

COMMON /COEFF/ COEFE( 8 , 10, 12 ) , COEFS( 8 , 6) ,NMBRSEG( 12) 
COMMON /COUNTER/ NODESAV, NODETOT, NBNODES, NPLANE 
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C 


COMMON /HEADER/ 
COMMON /INITA/ 
COMMON /INITB/ 
COMMON /INITC/ 
COMMON /NODNMBR/ 
COMMON /UNITS/ 


ITITLE( 20 ) , LINE 
I Z INDEX , MAPTEN , INCHES 

IBOX( 4,6), IRAY( 4,3) , NRAY, IPT1 ( 12) , IPT2( 12) 
PI ,RADDEG 

NODENUM( 4000) ,TOL(3) 

NU5,NU6,NU20 


DIMENSION NAME (20) 


C 

DATA PI /3. 141592654/ 

DATA RADDEG /57 . 29577951/ 

DATA IPT1 /l, 2, 4, 1,1, 2, 3, 4*5,6, 8, 5/ 

DATA IPT2 /2, 3, 3, 4, 5, 6, 7, 8, 6, 7, 7, 8/ 

DATA IRAY /I , 3 , 9 , 11 , 4 , 2 , 12 , 10 , 5 , 6 , 8 ,7/ 

DATA I BOX /4, 1,2, 3, 5, 1,6, 9, 12, 9, 10, 11, 8, 3, 7, 11, 5, 4, 8, 12, 6, 2, 7, 10/ 

C 



C INITIALIZE 



c 

NU5 * 5 
NU6 * 6 
NU20 - 20 
MUNIT » 20 
C 

NPLANE ■ 0 
NBNODES - 0 
NINODES * 0 
NODETOT - 0 
NODESAV - 0 


— COEFFICIENTS 


DO 20 1-1,8 

DO 10 J-1,6 
10 COEFS ( I , J ) - 0.0 

DO 20 J-1,10 
DO 20 K-1,12 
20 COEFE ( I , J , K) - 0.0 


NZONE - NUMBER OF 20NES TO BE USED TO GENERATE MESH 

I Z INDEX - NO OF THE ZONE TO BE COMPUTED 

- 1 DO ZONES 142 

- 3 DO ZONE 3 

- 4 DO ZONES 465 

- 6 DO ALL FIVE ZONES 
MAPTEN - 0 FIVE SEGMENTS PER EDGE 

« 1 TEN SEGMENTS PER EDGE 
INCHES « 0 COORDINATES INPUT IN FEET 

- 1 COORDINATES INPUT IN INCHES 


ECHO INPUT 

WRITE (NU6, 1100) 
C 

DO 90 1-1,500 
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READ ( NU5,1 000, END“ 100) NAME 
90 WRITE (NU6, 1110) NAME 
C 

100 REWIND NU5 
C 

C PRINT HEADER 

C 

WRITE (NU6, 11 20) 

C 

C READ PROBLEM DEFINITION 

C 

READ( NU5 , 1000 ) ITITLE 
READ(NU5 , 1010) N ZONE, I Z INDEX, 

& MAPTEN, INCHES 

C 

C WRITE PROBLEM DEFINITION 

C 

NRAY - 4 

C 

WRITE (NU 6, 1130) ITITLE 
WRITE (NU6 , 1140 ) N20NE , IZINDEX, 

& MAPTEN , INCHES 

C 

CALL RWIND( NU20 ) 

C 

C DRAW PICTURE 

C 

IDIM « 3 

CALL PICTURE(IDIM) 

C 

C HGM TWO DUCT DATA 

C 

CALL DCTDAT 

C 

RETURN 

C 

1000 FORMAT ( 20A4) 

1010 FORMATt 1615) 

1100 FORMATt 1H1,25X,17H INPUT DATA IMAGE // ) 

1110 FORMATt 5X, 2 0A4) 

1120 FORMATt 1H1 

3 / 40X, 34H LOCKHEED-HUNTSVILLE 

4 / 40X,34H VAX 11/785 VERSION 

5 / 40X,38H 2-DUCT HGM GEOMETRY MODULE ) 

1130 FORMATt/ 12H CASE TITLE: // 5X,20A4 ) 

1140 FORMATt// 21H GEOMETRY PARAMETERS: 

1 // 23H NZONE I 2 INDEX, 

2 22H MAPTEN INCHES ,// 6110 ) 

C 

END 

C 

c *. 

C ***** INPUT ************************************************* 
(;*»***»•»•••***** A************************ ****************** 

c 

SUBROUTINE INPUT 
C 

COMMON /COEFF/ COEFEt 8, 10, 12) ,COEFS( 8, 6) ,NMBRSEG( 12) 

COMMON /COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
COMMON /HEADER/ ITITLE ( 20) , LINE 
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C 


COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


/INITA/ I Z INDEX, MAPTEN, INCHES 

/INITB/ IBOX( 4,6) , I RAY ( 4,3) ,NRAY, IPT1( 12) , IPT2(12) 
/INPUTA/ EDGE( 3,12) ,POINT(3,8) ,SIDE(3,6) 

/MAPING/ MAPSIDE( 6 ) ,MAPSEG( 10,12) 

/MARCHS/ MARCH , INDEX ( 3 ) 

/MAXIMUM/ ETAMAX (10,12), SEGMAX (3,10,12) 

/NODNMBR/ NODENUM( 4000) ,TOL( 3) 

/SPACING/ ISTRTCH ( 3 ) , STRETCH ( 3 ) , ETAS (3,200), ETA{ 3 ) , DETA( 3 ) 
/UNITS/ NU5,NU6,NU20 

/ZONING/ I ZONE, ISECT,NZ INDEX, NMBRNDS( 3) 


DIMENSION IFLAGP( 8 ) , IFLAGE( 12 ) , IFLAGS( 6 ) 
DIMENSION LSTRTCH ( 3 ) ,LBCINPT(6) ,LMAPS(6) 
DIMENSION MAPEDGE (12), STRTCH ( 3 ) , LNMBRNDt 3 ) 


MAPEDGE(I) INDICATES TYPE OF GEOMETRY FOR EDGE I 

■ 1 LINEAR 

« 2 CIRCULAR ARC 

* 3 CONIC (PARABOLA, ELLIPSE, HYPERBOLA) 

- 4 HELICAL ARC 

« 5 TRIG FUNCTION OF X 

- 6 TRIG FUNCTION OF ANGLE 

■ 7 CUBIC POLYNOMIAL 


MAPS IDE ( I ) TYPE OF GEOMETRY FOR SURFACE I 

- 1 FLAT SURFACE 

- 2 CYLINDRICAL SURFACE 

- 4 EDGE OF REVOLUTION 

MARCH ETA DIRECTION IN WHICH COMPUTATION IS TO ADVANCE 


NMBRNDS ( I ) NUMBER OF NODES IN EACH ETA(I) DIRECTION 

ISTRTCH ( I ) « 0 NO STRETCHING IN ETA(I) DIRECTION 

■ 1 INPUT VALUES OF ETA ( I ) (N) 

- 2 DECREASE SPACING IN ETA(I) DIRECTION (INPUT FACTOR) 

- 3 INCREASE SPACING IN ETA(I) DIRECTION (INPUT FACTOR) 

■ 4 DOUBLE STRETCHING (INPUT FACTOR) 

- 5 DECREASING SPACING IN ETA(I) DIRECTION (MINIMUM) 

- 6 INCREASING SPACING IN ETA(I) DIRECTION (MINIMUM) 

■ 7 DOUBLE STRETCHING (MINIMUM SPACING) 

■ 8 ORIGINAL DECREASING STRETCHING FUNCTION 

- 9 ORIGINAL INCREASING STRETCHING FUNCTION 
-10 USER INPUT STRETCHING FUNCTION 


WRITE (NU6, 1100) ITITLE , IZONE 


C 

C 

C 

c 

c 


READ INPUT PARAMETERS 
100 IF(MAPTEN.EQ.O) THEN 

READ(NU5, 1000) 
END IF 

READ(NU5, 1000) 
READ( NU5, 1000 ) 

IF ( MARCH . LT. 1 .OR. MARCH. GT. 3) 


MAPEDGE 

MAPS IDE 
MARCH 

MARCH - 1 
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READ (NU 5, 1000) NMBRNDS , ISTRTCH 


C 

C PRINT OUT INPUT VARIABLES 

C 

WRITE (NU6, 3000) MAPEDGE 
WRITE (NU6, 3010) MAPSIDE 
WRITE (NU6, 3020) MARCH 

C 

WRITE (NU6, 3040) NMBRNDS , ISTRTCH 

C 

C INITIALIZE FLAGS TO CHECK IF POINT, EDGE, OR SURFACE HAS BEEN INPUT 

C 

DO 130 1-1,8 
130 IFLAGP(I) - 0 

C 

DO 140 1-1,12 
140 IFLAGE(I) - 0 


DO 150 1-1,6 
150 I FLAGS ( I ) - 0 
C 

C SET ORDER OF EXECUTION 

C 

INDEX ( 1 ) > MARCH 
INDEX( 3 ) - 3 
C 

DO 160 1-2,3 

INDEX(I) • INDEX(I-l) + 1 
160 I F ( INDEX ( I) .GT.3) INDEX(I) - 1 
C 

C INITIALIZE ETA 

C 

ETA(l) - 0.0 
ETA( 2 ) - 0.0 
ETA( 3 ) - 0.0 
C 

C — — 

C READ STRETCHING PARAMETERS 

C 

c 

300 DO 310 1-1,3 
C 

STRETCH ( I ) - 0.0 
C 

I F ( ISTRTCHt I ) . LE. 1 ) GO TO 310 

C 

READ(NU5, 1020) STRETCH 

C 

GO TO 320 
C 

310 CONTINUE 



C COMPUTE STRETCHING FUNCTION PARAMETER B USING NEWTON-RAPHSON 



320 DO 370 1-1,3 
C 

B - STRETCH ( I ) 

DS - B 

TNODE - REAL ( NMBRNDS ( I ) ) 

C 
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DECREASING OR INCREASING SPACING (INPUT MINIMUM SPACING) 

IF( ISTRTCH( I) . EQ. 5 .OR. ISTRTCH ( I ) . EQ . 6 ) THEN 

B » SORT (1.0 - ( TNODE - 1.0)*DS)/(TNODE - 1.0) 

DO 330 ITER»1,10 

ARG1 ■ B* (TNODE - 1.) 

ARG2 » B*( TNODE - 2.) 

EXPAl - EXP(ARGl) 

EXPA2 • EXP( ARG2) 

TANHl » (EXPAl - 1. /EXPAl )/( EXPAl + 1. /EXPAl) 

TANH2 ■ (EXPA2 - 1 ./EXPA2)/(EXPA2 ♦ 1./EXPA2) 

PSI - DS - (1.0 - TANH2/TANH1) 

TANHPl ■ 1.0 - TANHl** 2 
TANHP2 « 1.0 - TANH2**2 

PSIP « ( 1 . 0/TANH1 ) * ( TANHP 2* ( TNODE - 2.) 

6 - ( TANH2/TANH1 ) * TANHPl* (TNODE - 1.)) 

IF( PSIP. EQ. 0 . ) THEN 

WRITE (NU6, 5100) 

STOP 

END IF 


B0 » B 

B « B0 - PSI/PSIP 
DBF » (B - B0)/B0 

WRITE (NU6, 5200) I , PSI , PSIP, B, DBF 

330 IF(ABS(DBF) .LE. 0.001) GO TO 340 

340 END IF 

DOUBLE STRETCHING (INPUT MINIMUM SPACING) 


C 

C 

c 


c 


c 


c 


IF(ISTRTCH(I).EQ.7) then 

B - SORT (1.0 - (TNODE - 1.0)*DS)/(TNODE - 1.0) 

DO 350 ITER-1, 10 

ARG1 - B*( TNODE >1.) 

ARG3 - B* (TNODE - 3.) 

EXPAl • EXP(ARGl) 

EXPA3 - EXP( ARG3 ) 

TANHl « (EXPAl - 1. /EXPAl )/( EXPAl + 1. /EXPAl) 

TANH3 » (EXPA3 - l./EXPA3)/(EXPA3 + 1./EXPA3) 
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C 


C 

& 

c 

c 

c 

c 

c 


c 

c 

350 

c 

360 

C 

370 STRETCH ( I ) - B 


PSI - DS - 0.5* (1.0 - TANH3/TANH1 ) 

TANHPl ■ 1.0 - TANHl**2 
TANHP3 » 1.0 - TANH3**2 

PSIP » ( 0 . 5/TANH 1 ) * ( TANHP 3 * ( TNODE - 3.0) 

- (TANH3/TANH1 ) * TANHPl* (TNODE - 1.0)) 

IF(PSIP.EQ.O. ) THEN 

WRITE(NU6,5000) 

STOP 

END IF 


B0 « B 

B » B0 - PSI/PSIP 
DBF - (B - B0)/B0 

WRITE (NU6, 5200) I, PSI, PSIP, B, DBF 

IF(ABS(DBF) .LE. 0.001) GO TO 360 

END IF 


—WRITE STRETCHING PARAMETERS 
DO 390 1-1,3 

IF ( ISTRTCH( I ) . LE . 1 ) GO TO 390 

WRITE(NU6, 3050) STRETCH 
GO TO 400 

390 CONTINUE 

READ INPUT PARAMETERS FOR ARBITRARY GRID SPACING (ETAS) 
400 DO 410 1-1,3 

IF(ISTRTCH(I).NE.l) GO TO 410 

IF( ISECT.GT.l .AND. I. NE. MARCH) GO TO 410 

READ(NU5, 1020) ( ETAS ( I , J ) , J-l , NMBRNDS( I ) ) 

WRITE ( NU6 , 1110 ) I , ( ETAS ( I , J ) , J- 1 , NMBRNDS ( I ) ) 

410 CONTINUE 

HGM TWO DUCT INPUT DATA 

CALL HGMIN 

INPUT PARAMETERS FOR EDGE COEFFICIENTS 

WRITE (NU6, 1100) ITITLE,IZONE 
WRITE ( NU 6, 2040) 
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LINE 

II 


6 

0 


TOTAL NUMBER OF EDGES 

NEDGES ■ 8*3 - 12 

DO 540 1-1, NEDGES 

IF(IFLAGE(I) .EQ.l) GO TO 540 

ITOTAL « 1 

MAP - MAPEDGE ( I ) 

DETERMINE THE NUMBER OF SEGMENTS ON AN EDGE 

DO 500 J«l,10 


NMBRSEG(I) 

MAP 


J 

MAP/ 10 


I F ( MAP . EQ • 0 ) GO TO 510 

500 ITOTAL * ITOTAL*10 

510 MAP » MAPEDGE ( I ) 

DO 530 J»1,NMBRSEG(I) 

DETERMINE THE MAPPING FOR EACH SEGMENT 

MAPSEG( J , I ) ■ MAP/ITOTAL 

MAP - MAP - MAPSEG(J,I)*ITOTAL 

EDGE COEFFICIENTS FOR EACH SEGMENT 

I F ( MAPSEG (J.I).LE.l) GO TO 530 

IF(I.NE.II) THEN 

LINE * LINE + 2 

IFUINE.GE.60) THEN 


C 

c 

c 

c 

c 

c 

c 


WRITE (NU6, 11 00) ITITLE.IZONE 
WRITE (NU6, 2040) 

LINE - 8 

END IF 

WRITE(NU6,2050) I,J, (COEFE(K, J, I) , K-l ,8 ) 
ELSE 

LINE - LINE + 2 
IFUINE.GE.60) THEN 
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C 

C 

C 

C 

C 

C 

C 


WRITE (NU6, 1100) ITITLE,IZONE 
WRITE (NU6, 2040) 

LINE * 8 

END IF 

WRITE(NU6, 2060) J, (C0EFE(K,J,I) ,K»1,8) 

END IF 

II - I 

530 ITOTAL ■ ITOTAL/10 
540 CONTINUE 

INPUT PARAMETERS FOR SURFACE COEFFFICIENTS 


LINE ■ LINE + 5 
IF( LINE .GE . 60 ) THEN 

WRITE (NU6, 11 00) ITITLE,IZONE 
LINE « 6 
END IF 


WRITE (NU6, 3060) 


C 

C 


DO 610 1-1,6 

I F ( I FLAGS (I).EQ.l) GO TO 610 

IF(MAPSIDE(I) .LE.2) GO TO 610 

LINE - LINE + 2 

IF( LINE.GE.60) THEN 

WRITE (NU6, 1100) ITITLE,IZONE 
WRITE (NU6, 3060) 

LINE - 8 

END IF 

WRITE (NU6, 3 06 5) I , (COEFS( J , I) , J-l , 8) 

610 CONTINUE 

INPUT DATA FOR CORNER POINTS AND SEGMENT END POINTS 
LINE - LINE + 5 
IFUINE.GE.60) THEN 

WRITE (NU6, 11 00) ITITLE,IZONE 
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LINE » 6 
END IF 

WRITE (NU6, 3070) 

DO 750 I«1,NEDGES 
I F ( I . GT . 8 ) GO TO 700 

CORNER POINTS 

LINE - LINE + 2 
IF ( LINE . GE . 60 ) THEN 

WRITE (NU6, 1100) ITITLE,IZONE 
WRITE (NU6, 3070) 


C 

C 

C 

C 

C 


LINE ■ 8 
END IF 

WRITE (NU6, 3080) I , ( POINT{ J , I) , J«1 , 3) 


700 IF(NMBRSEG(I) .EQ.l) GO TO 750 
C 

C SEGMENT MAXIMUMS 

C 

DO 740 J-1,NMBRSEG(I) - 1 

C 

LINE « LINE + 2 

C 

IF( LINE. GE. 60) THEN 
C 


WRITE (NU6# 1100) ITITLE.IZONE 
WRITE (NU6, 3070) 


C 

C 

C 

6 

C 

c 


c 

720 

C 

730 

C 

C CONVERT 

C 

C 


LINE - 8 
END IF 

IF(IFLAGE(I) .EQ.l) ETAMAX(J,I) * 

ETAMAX ( J , I ) * ( NMBRNDS ( N ) - 1.0) + 1.0 

WRITE ( NU6 , 3090 ) I » J » ( SEGMAX(K, J , I) ,K»1, 5) ,ETAMAX( J, I) 

DO 720 N-1,3 
DO 720 L-1,4 

IF(IRAY(L,N) .EQ.l) GO TO 730 
CONTINUE 

NODE NUMBER TO ETA VALUE 

ETAMAX(J,I) ■ (ETAMAX(J,I) - 1.0)/(NMBRNDS(N) - 1.0) 
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740 CONTINUE 
750 CONTINUE 
C 

DO 930 1-1,3 
C 

DETA(I) - 1 . 0/ ( NMBRNDS ( I ) - 1.0) 

C 

FACTOR - NMBRNDStnnO.O 
IF(ISTRTCH(I) .GT.0) FACTOR - NMBRNDS ( I ) * 20 . 0 

C 

DO 920 J-l , NRAY 

C 

L - IRAY( J , I) 

— COMPARE X, Y, AND Z BETWEEN END POINTS OF AN EDGE 
DO 920 K-1,3 

DELTA - ABS(POINT(K, IPTl(L) ) - POINT( K, IPT2(L) ) )/FACTOR 

IF ( DELTA. LE .0.0) GO TO 920 

IF( DELTA. LT. TOL( K) ) TOL(K) « DELTA 

920 CONTINUE 
930 CONTINUE 


PRINT TITLE 


WRITE (NU 6, 1100) ITITLE,IZONE 
LINE - 1 
WRITE (NU6, 1140) 
LINE - 3 


RETURN 

FORMAT STATEMENTS 

1000 FORMAT ( 1615) 

1010 FORMAT ( 6110 ) 

1020 FORMAT ( 8E10 . 0 ) 

1030 FORMAT( (4110) ) 

C 

1100 FORMAT( 1H1 , 10X, 20A4 , 13X, 6H ZONE, 13) 

1110 FORMAT ( // 11H FIXED ETA ,I1,8H VALUES: // ( 10 { 3X,F10.7) ) ) 

1140 FORMAT ( / 44H NODE X Y Z ) 

C 

2000 FORMAT (// 32H EDGE SHAPE FUNCTION INDICATORS: 

6 // 40H EDGE 1 EDGE__2 EDGE_3 EDGE_4 / 4110 ) 

2010 FORMAT (// 26H BOUNDARY CONDITION FLAGS: 

1 // 40H EDGE 1 EDGE 2 EDGE 3 EDGE_4 / 4110 

2 /// 26H MARCHINS DIRECTION - ETA ,lT) 

2020 FORMAT ( // 17H NUMBER OF NODES: 

1 // 20H ETA 1 ETA 2 / 2110 

2 // 22H STRETCHING FUNCTIONS: 

3 // 20H ETA 1 ETA 2 / 2110 ) 

2030 FORMAT ( // 32H STRETCHING FUNCTION PARAMETERS: 


C-ll 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 


ooooo oonno 


LMSC-HEC TR F225968 


* // 

2040 FORMAT (// 

1 // 

2 
3 

2050 FORMAT ( / 
2060 FORMAT (/ 
2080 FORMAT*/ 
2090 FORMAT*/ 

C 


30H ETA_1 ETA_2 / 2(2X,F13.7) ) 

25H EDGE SHAPE COEFFICIENTS! 

42H EDGE SEGMENT COEFF 1 COEFF 2, 

42H COEFF_3 COEFF“4 COEFF”5, 

42H COEFF_6 COEFF~7 COEFF~8 ) 

2X,I2,6X,I1,3X,8(1X,F13.7) ) 

10X,I1,3X,8(1X,F13.7)) 

7H POINT ,I2,1H: ,12X,2(7X,F13.7) ,9X,F7.2) 

7H EDGE ,I2,1H: , 7X, 13 , 2X, 2 ( 7X, F13 . 7 ) , 9X,F7 . 2 , 9X, F9 . 2 ) 


3000 FORMAT*// 
1 // 
2 
3 

3010 FORMAT*// 
1 // 
2 

3 / 

3020 FORMAT*// 
3040 FORMAT*// 
1 // 
2 // 
3 // 

3050 FORMAT*// 

1 // 

2 / 
3060 FORMAT*// 

1 // 
2 
3 

3065 FORMAT* / 
3070 FORMAT*// 
1 // 
2 

3080 FORMAT*/ 
3090 FORMAT*/ 

& 


ETA_3 / 3110 


32H EDGE SHAPE FUNCTION INDICATORS: 

40H EDGE_1 EDGE_2 EDGE_3 EDGE_4 , 

40H EDGE_5 EDGE 6 EDGE 7 EDGE 8, 

4 OH EDGE_9 EDGE_l0 EDGE_Il EDGE_l2 / 12110 ) 

35H SURFACE SHAPE FUNCTION INDICATORS: 

45H SURFACE_1 SURFACE_2 SURFACE 3, 

45H SURFACE 4 SURFACE_5 SURFACE~6 

6(10X,I?)) 

26H MARCHING DIRECTION ■ ETA_,Il) 

17H NUMBER OF NODES: 

3 OH ETA_1 ETA_2 

22H STRETCHING FUNCTIONS: 

3 OH ETA_1 ETA_2 ETA 3 / 3110 ) 

32H STRETCHING FUNCTION PARAMETERS: 

45H ETA 1 ETA 2 

3(2X,F1377) ) 

28H SURFACE SHAPE COEFFICIENTS: 

42H SURFACE COEFF 1 

42H COEFF 3 COEFF~4 

42H COEFF 6 COEFF“7 

4X,Il,9X,8(lX,Fl3.7) ) 

29H COORDINATES AND ETAMAXES: 

10X, 39H SEGMENT X Y, 

41H Z ETAMAX ) 

7H POINT ,12, 1H: , 8X, 3 ( 6X, F13 . 7) , 2* 8X, F7 . 2 ) ) 

7H EDGE ,I2,1H: , 3X, 13 , 2X, 3 ( 6X, F13 . 7 ) ,2(8X,F7.2) 

,7X,F9. 2) 


ETA 3 


COEFF 2, 
COEFF“5, 
COEFF”8 ) 


C 

5000 FORMAT* 1H0, 4 3H SINGULARITY IN TWO-END STRETCHING FUNCTION, 
& 10H PARAMETER) 

5100 FORMAT* 1 HO, 4 3H SINGULARITY IN ONE-END STRETCHING FUNCTION, 
& 10H PARAMETER) 

5200 FORMAT(1HO,4H I -,I2,7H PSI «,E12.4,8H PSIP -,E12.4, 

6 5H B -,E12.4,7H DBF -,E12.4) 

C 

END 


*****MAPPING 


SUBROUTINE MESH 


MESH CONTROLS THE GENERATION OF THE SPATIAL COORDINATES OF EACH NODE 
IN THE MESH. 


COMMON /COUNTER/ NODESAV, NODETOT , NBNODES , NPLANE 
COMMON /HEADER/ ITITLE ( 20 ) , LINE 
COMMON /INITA/ IZINDEX, MAPTEN, INCHES 
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COMMON /INITB/ 
COMMON /INPUTA/ 
COMMON /INPUTBC/ 
COMMON /MAPING/ 
COMMON /MARCHS/ 
COMMON /NODNMBR/ 
COMMON /OUT/ 
COMMON /PARTIAL/ 
COMMON /SPACING/ 
COMMON /UNITS/ 
COMMON /ZONING/ 


IBOX(4,6) ,IRAY(4,3) ,NRAY, IPT1 ( 12) , IPT2 ( 12 ) 

EDGE( 3,12) r POINT ( 3 ,8), SIDE (3, 6) 

INODEBC( 3 ) , ISIDE( 3 ) 

MAPS IDE ( 6 ) , MAPS EG ( 10,12) 

MARCH , INDEX( 3 ) 

NODENUM( 4000 ) , TOL( 3 ) 

NODE (3, 4000) 

DEDN (3,12), DSDN (3,2), SNORMAL (3,6) 

ISTRTCH( 3) , STRETCH ( 3 ) , ETAS (3,200) »ETA( 3) , DETA( 3) 
NU5,NU6,NU20 

I ZONE , ISECT , NZ INDEX , NMBRNDS ( 3 ) 


DIMENSION E( 26) 

DIMENSION IFACE1 (3,3), IFACE2( 3,3), IS I DBS ( 6 ) 
DIMENSION TANGENT ( 3 ) ,VECTOR( 3) ,DIRECT(3) 


C 

DATA IFACE1 /0 , 1 , 2 , 1 , 0 , 5 , 2, 5 , 0/ 
DATA IFACE2 /0 , 3 , 4 , 3 , 0 , 6 , 4 , 6 , 0/ 
DATA ISIDES /3 , 4 , 1 , 2 , 6, 5/ 

C 

INIT - 1 


C 


FOLD « 0 


C 

C DETERMINE SURFACES 

C 

ISIDE1 * IFACE1 (MARCH* INDEX ( 2) ) 
ISIDE2 * IFACE2( MARCH, INDEX ( 2) ) 
ISIDE3 - IFACE1 ( MARCH , INDEX ( 3 ) ) 
ISIDE4 » I FACE 2 ( MARCH , INDEX ( 3 ) ) 
ISIDE5 - IFACEl ( INDEX( 2 ) , INDEX ( 3 ) ) 
ISIDE6 - IFACE2( INDEX( 2 ) , INDEX ( 3 ) ) 


C 

C— 

C 


C- 

C 

C- 


-MAXIMUM NUMBER OF NODES IN A PLANE 

MAXPL - NMBRNDS ( INDEX( 2) ) ‘NMBRNDS ( INDEX ( 3) ) 

INDEX ( I ) - MARCH 
AXIS 

DO 700 IAXIS*1 , NMBRNDS ( INDEX ( 1 ) ) 

IF( IAXIS.EQ.l .AND. ISECT. GT.l) GO TO 700 
— SEPERATE BOUNDARY CONDITIONS AND DETERMINE ETA 
CALL ETABC { MARCH, INDEX (1),IAXIS) 

—CALCULATE COORDINATES AND DERIVATIVES FOR POINTS ON EDGES 
CALL EDGES ( INIT , INDEX( 1 ) , ETA ( INDEX ( 1 ) ) , FOLD ) 

—NUMBER OF NODES TO BE STORED 
NODSTOR - NODESAV 
_ 


c 

c- 

c 

c 

c- 

c 

c 

c- 

c 

c 

c- 


C ROW 
c 
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DO 500 JAXIS-1,NMBRNDS(INDEX(2)) 

— SEPERATE BOUNDARY CONDITIONS AND DETERMINE ETA 
CALL ETABC( MARCH* INDEX( 2) , JAXIS) 

— CALCULATE COORDINATES AND DERIVATIVES FOR POINTS ON EDGES 
CALL EDGES ( INIT, INDEX( 2) , ETA( INDEX( 2) ) * FOLD) 

— CALCULATE COORDINATES AND SURFACE NORMAL FOR POINTS ON SURFACES 
CALL SURFACE ( INIT, ISIDE1) 

CALL SURFACE (INIT, IS IDE 2 ) 


COLUMN 


DO 400 KAXIS“1 , NMBRNDS ( INDEX ( 3 ) ) 

—SEPERATE BOUNDARY CONDITIONS AND DETERMINE ETA 
CALL ETABC (MARCH, INDEX( 3 ) , KAXIS ) 

— CALCULATE COORDINATES AND DERIVATIVES FOR POINTS ON EDGES 


CALL EDGES ( INIT , INDEX ( 3 ) , ETA( INDEX( 3 ) ) , FOLD) 


CALCULATE COORDINATES AND SURFACE NORMAL FOR POINTS ON SURFACES 

CALL SURFACE ( INIT , ISIDE3 ) 

CALL SURFACE ( INIT, ISIDE4 ) 

CALL SURFACE ( INIT, ISIDE5) 

CALL SURFACE( INIT, ISIDE6) 


ETA COEFFICIENTS FOR TRI-LINEAR INTERPOLATION 


E(l) 

m 

1.0 - ETA( 3) 

E ( 2 ) 

m 

ETA ( 3 ) 

E(3) 

m 

1.0 - ETA( 2 ) 

E( 4 ) 

u 

ETA( 2) 

E( 5) 

m 

1.0 - ETA(l) 

E( 6 ) 

m 

ETA( 1 ) 

E(7) 

m 

E(5)*E(3) 

E( 8 ) 

m 

E( 5 ) *ETA( 2) 

E( 9 ) 

m 

ETA( 1 ) *E( 3 ) 

E( 10 ) 

m 

ETA( 1 ) *ETA( 2 ) 

E( 11) 

m 

E( 5 ) *E( 1 ) 

E( 12 ) 

m 

E( 5) *ETA( 3) 

EC 13) 

m 

ETA( 1 ) *E(1 ) 

E{ 14 ) 

m 

ETA( 1 ) *ETA( 3) 

E( 15) 

m 

E(3)*E(1) 

E( 16) 

m 

E(3)*ETA(3) 

E( 17) 

• 

ETA(2)*E(1) 

E(18) 

m 

ETA(2)*ETA(3) 

E( 19 ) 

m 

E(5)*E(3) 

E{ 20) 

* 

E(5)*E(3) 


*E(1) 
*ETA( 3) 
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E( 21 ) » E( 5) *ETA( 2) *E( 1 ) 

E( 22) - E( 5) *ETA( 2 ) *ETA( 3) 
E( 23 ) - ETA( 1 ) *E( 3 ) *E(1) 

E ( 2 4 ) - ETA( 1 ) *E( 3 ) *ETA(3) 

E( 25 ) - ETA ( 1 ) * ETA ( 2 ) * E ( 1 ) 

E( 26) ■ ETA ( 1 ) * ETA { 2 ) * ETA ( 3 ) 


C INCREMENT NODE COUNTERS 


C 


INIT - 0 


C 

NODESAV » NODESAV + 1 
NODNUM » NODESAV + NODETOT 
C 


C CALCULATE THE COORDINATES 

C 


DO 340 L«l,3 
340 NODE (L, NODESAV) ■ 

1 + 

2 + 

3- E ( 7 ) *EDGE ( L, 5 ) -E( 

4- E(ll)*EDGE( 

5- E(15)*EDGE( 

6 

7 + 

8 + 

9 + 

350 CONTINUE 


E(1)*SIDE(L,1) + 
E( 3 ) *SIDE( L f 2 ) + 
E(5)*SIDE(L,5) + 
8 ) *EDGE( L, 8 )-E( 


E( 2 ) *SIDE( L, 3 ) 

E( 4 ) *SIDE( L, 4 ) 

E( 6 ) *SIDE( L, 6 ) 

9 ) *EDGE( L, 6 )-E( 10 ) *EDGE ( L, 


L,5)-E( 8 ) *EDGE( L, 8)-E( 9 ) *EDGE( L, 6 )-E( 10 ) ‘EDGE ( L, 7) 
L, 4 )-E( 12 )*EDGE( L, 12)-E( 13)*EDGE( L, 2)-E( 14 ) *EDGE( L, 10 ) 
L, 1 )-E( 16 ) *EDGE( L, 9 ) -E( 17 ) *BDGE(L, 3 ) -E( 18 ) *EDGE( L, 11 ) 
+ E( 19 )* POINT ( L, 1 ) + E(20)* POINT (L f 5) 

n / \ T Mm / r A \ x B \ * AA T UW t V O \ 


E(21)*POINT(L,4) + 
E(23)*POINT(L,2) + 
E( 25) * POINT (L, 3 ) + 


E( 22 ) *POINT( L,8 ) 
E( 24)*POINT(L,6) 
E(26)*POINT(L,7) 


400 CONTINUE 
500 CONTINUE 


C 

C 

c 

c 

c 

c 

c 

c 


PLANE OUTPUT 

IF ( IAXIS . EQ. 1 ) GO TO 700 

-PRINT AND STORE DATA 

IF ( NODSTOR. NE • 0 ) THEN 

CALL OUTPUT (NU 20, NODSTOR) 


C TRANSFER DATA SECOND PLANE TO FIRST PLANE 


C 

C 

600 

C 

630 

C 


DO 630 I»l, NODESAV 
DO 600 J-1,5 

NODE ( J r I ) * NODE( J , I + NODSTOR) 
CONTINUE 


END IF 

C 

C TRANSFER NODE NUMBERS FROM SECOND PLANE TO FIRST PLANE 

C 

DO 650 L»1,MAXPL 

650 NODENUM(L) - NODENUMU + MAXPL) 


C 

700 CONTINUE 
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: 

OUTPUT 






C— — 

c 

c— 

-PRINT AND 

STORE DATA 








NODSTOR « 

NODESAV 


C 

I F ( NODES AV . NE . 0 ) CALL OUTPUT(NU20, 

NODSTOR) 


c 




NODESAV - 

0 



L 

p 

RETURN 







L 

c— 

-FORMAT STATEMENTS 





1100 FORMAT ( 1H1 , 10X, 20A4 , 13X, 8H SECTION 

,I2,3H OF 

,I3,9H FOR ZONE, 13) 

V. 

1110 FORMAT( 


8H NODE 

«,I10 





1 

/ 

41X,2H X, 

13X,2H Y, 12X 

x 

CM 

Z 



2 

/ 

8H E( 1) 

»,F13.7, 15H 


SIDE(l) 

»,6(1X,F13.7) 


3 

/ 

8H E ( 3) 

-,F13.7,15H 


SIDE( 2 ) 

-,6(1X,F13.7) 


4 

/ 

8H E ( 2) 

■,F13.7,15H 


SIDE( 3 ) 

-,6(1X,F13.7) 


5 

/ 

8H E( 4) 

■ ,F13 . 7, 15H 


SIDE( 4 ) 

■,6(1X,F13.7) 


6 

/ 

8H E ( 5) 

«,F13.7,15H 


SIDE( 5) 

»,6(1X,F13.7) 

n 

7 

/ 

8H E ( 6) 

■ f F13.7,15H 


SIDE(6) 

»,6(1X,F13.7) ) 

1120 FORMAT ( 


8H E( 15 ) 

-,F13.7,15H 


EDGE( 1) 

-,6(1X,F13.7) 


1 

/ 

8H E( 13) 

»,F13.7,15H 


EDGE( 2) 

»,6( IX, F13.7) 


2 

/ 

8H E ( 17 ) 

«,F13.7,15H 


EDGE( 3) 

-,6(1X,F13.7) 


3 

/ 

8H E ( 11 ) 

-,F13.7,15H 


EDGE( 4) 

«,6(1X,F13.7) 


4 

/ 

8H E( 7) 

■»F13.7,15H 


EDGE( 5) 

*,6(1X,F13.7) 

/■» 

5 

/ 

8H E ( 9) 

-,F13.7,15H 


EDGE( 6) 

■ , 6( 1X,F13 .7 ) ) 

1130 FORMAT ( 

8H E ( 10 ) 

-,F13.7,15H 


EDGE( 7) 

■ , 6( IX, F13. 7 ) 


1 

/ 

8H E ( 8) 

«,F13.7,15H 


EDGE( 8) 

«,6(1X,F13.7) 


2 

/ 

8H E( 16) 

-,F13.7,15H 


EDGE( 9) 

»,6(1X,F13.7) 


3 

/ 

8H E ( 14 ) 

■ ,F13 .7 , 15H 


EDGE ( 10) 

-,6(1X,F13.7) 


4 

/ 

8H E( 18 ) 

■ r F13 . 7 , 15H 


EDGE(ll) 

»,6(1X,F13.7) 

n 

5 

/ 

8H E ( 12 ) 

-,F13.7,15H 


EDGE ( 12) 

-,6(1X,F13.7)) 

1140 FORMAT ( 

8H E ( 19 ) 

■»F13.7, 15H 


POINT ( 1 ) 

-,6(1X,F13.7) 


1 

/ 

8H E ( 23 ) 

« ^ F13 • 7 , 15H 


POINT ( 2 ) 

■ , 6 ( IX, F13 .7 ) 


2 

/ 

8H E( 25) 

-,F13.7,15H 


POINT ( 3 ) 

-,6(1X,F13.7) 


3 

/ 

8H E( 21) 

■,F13.7,15H 


POINT( 4) 

«,6(1X,F13.7) 


4 

/ 

8H E( 20) 

■»F13.7,15H 


POINT ( 5 ) 

«,6(1X,F13.7) 


5 

/ 

8H E ( 24 ) 

■#F13 .7 , 15H 


POINT ( 6 ) 

»,6( 1X,F13.7) 


6 

/ 

8H E ( 26 ) 

»,F13.7,15H 


POINT(7) 

-,6(1X,F13.7) 

P 

7 

/ 

8H E( 22) 

«,F13.7,15H 


POINT ( 8 ) 

»,6(1X,F13.7) ) 

c 

Q*i 

END 

»* 



***) 



C*****MAPPING* 



*** 



C 

SUBROUTINE EDGES ( INIT, IDIR, ETA 

, FOLD) 


c— 

c 

INTERPOLATES ALONG 

EACH EDGE OF 

EACH SIDE OF 

ZONE AS PART OF THE 

c 

c— 

BI/TRI-LINEAR INTERPOLATION SCHEME 
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COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


/COEFF/ COEFE(8,10,12) ,COEFS( 8 , 6) ,NMBRSEG( 12 ) 

/COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
/HEADER/ ITITLE(20),LINE 
/INITA/ I Z INDEX, MAPTEN, INCHES 

/INITB/ I BOX ( 4 , 6 ) , I RAY ( 4 ,3) , NRAY, IPT1 ( 12) , IPT2( 12 ) 

/INPUTA/ EDGE ( 3 ,12) ,POINT( 3 ,8), SIDE (3, 6) 

/MAPED/ KSEG( 12 ) ,UAXIS( 3,6), I EDGE 1 (6) , IEDGE2 ( 6 ) , GAMMA 

/MAXIMUM/ ETAMAX( 10,12) , SEGMAX( 3,10,12) 

/PARTIAL/ DEDN (3,12) , DSDN (3,2) ,SNORMAL( 3 , 6 ) 

/SEGMENT/ PT1 ( 6 , 10, 12) ,PT2 (6,10,12), ETA1 (10,12), ETA2( 10,12) 
/UNITS/ NU5,NU6,NU20 

/ZONING/ I ZONE , ISECT , NZ INDEX , NMBRNDS ( 3 ) 


INTERMEDIATE PRINT 


NODNUM » NODESAV + NODETOT + 1 


CALCULATE EDGE COORDINATES AND DERIVATIVE 
20 DO 200 1*1, NRAY 
— DETERMINE WHICH EDGE 

IEDGE * IRAY( I , IDIR) 

— EDGE INITIALIZATION 

IF(INIT.EQ.l) CALL EMAPI ( IEDGE , FOLD, IDIR, I ) 

— DETERMINE WHICH SEGMENT 

DO 100 JSEG-1,NMBRSEG( IEDGE) 

ISEG -JSEG 

100 IF ( ETA. GE.ETAM JSEG, IEDGE) .AND. ETA. LE. ETA2( JSEG, IEDGE) ) GOTO 110 

110 KSEG( IEDGE) * ISEG 

DETERMINE WHERE ALONG THE SEGMENT 

DENOM - ETA2( ISEG, IEDGE) - ETAl ( ISEG, IEDGE) 

IF ( DENOM . EQ . 0 . 0 ) DENOM - 1.0 

RATIO - (ETA - ETA1( ISEG, IEDGE) ) /DENOM 

— CALCULATE THE COORDINATES AND DERIVATIVE 

CALL EMAP(IEDGE, ISEG, RATIO, EDGE(1, IEDGE) ,DEDN(1, IEDGE) ) 


INTERMEDIATE PRINT 


200 CONTINUE 
RETURN 

— FORMAT STATEMENTS 
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C 


1100 FORMAT(1H1,10X,20A4,13X,8H SECTION, 12 , 3H OF,I3,9H FOR ZONE, 13) 
1110 FORMAT(/ 41H NODE EDGE X Y, 

1 43H Z TANGENTS THETA, 

2 36H PHI ETA RATIO ETAMAX) 

1120 FORMAT ( 1X,I6,4X,I2,3X,3(3X,F13.7),11X,2(2X,F7.2),3(2X,F7.5)) 


END 


****MAPPING 


SUBROUTINE EMAPI ( IEDGE, FOLD, IDIR,NMBREDG) 
EDGE MAPPING INITIALIZATION 


COMMON 

COMMON 

COMMON 

& 

COMMON 

& 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


/COEFF/ COEFE (8,10,12), COEFS (8,6), NMBRSEG ( 12 ) 

/COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
/EDGE0/ UI(3,10,12),UJ(3,10,12),UK(3,10,12) , 
R1(3,10,12),R2(3,10,12) , THETA ( 10,12) 

/EDGE3/ ARC (10,12), ARC 1(10,12) ,XLENGTH( 10,12) , 

RA ( 1 0 , 1 2 ) , RC ( 1 0 , 1 2 ) , RE (10, 12) ,THETAl (10,12) 
/EDGE8/ RM1 (10,12) ,RM2 (10,12) ,RPAl(10,12) ,RPA2(10,12) 
/HEADER/ ITITLE ( 20 ) , LINE 
/INITA/ I Z INDEX, MAPTEN, INCHES 

/INITB/ IBOX( 4,6), IRAY( 4 , 3 ) ,NRAY, IPT1 ( 12) , IPT2( 12) 

/INITC/ PI ,RADDEG 

/INPUTA/ EDGE (3, 12) , POINT (3,8) , SIDE (3,6) 

/MAPED/ KSEG ( 12 ) , UAXIS (3,6), IEDGEl ( 6 ) , IEDGE2 { 6 ) , GAMMA 

/MAPING/ MAPSIDE(6) ,MAPSEG( 10,12) 

/MARCHS/ MARCH , INDEX ( 3 ) 

/MAXIMUM/ ETAMAX (10,12) ,SEGMAX{ 3,10, 12) 

/PARTIAL/ DEDN(3,12) ,DSDN( 3,2) ,SNORMAL(3,6) 

/SEGMENT/ PT1 (6,10,12) ,PT2( 6, 10 , 12) ,ETAl ( 10, 12) ,ETA2( 10, 12 ) 


DIMENSION ETAMX (10,12) ,VECTER(3) 


— INITIALIZE 

KSEG( IEDGE) - 1 

DO 10 N-l, NMBRSEG (IEDGE) 

10 ETAMX (N, IEDGE) ■ ETAMAX ( N, IEDGE) 

40 ETAMX ( NMBRSEG ( IEDGE) , IEDGE) • 1.0 

ETAM1, IEDGE) ■ 0.0 
ETAMAX ( NMBRSEG ( IEDGE) , IEDGE) ■ 1.0 


INITIALIZE SEGMENTS 

DO 1000 ISEG-1, NMBRSEG ( IEDGE) 

DETERMINE THE COORDINATES AT THE END POINTS 

IF( ISEG.EQ.l) THEN 
C 

DO 50 J-1,3 


PT1(J, 1, IEDGE) ■ POINT (J,IPTl( IEDGE)) 
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50 

C 

c 


60 

c 

c 

c ETA value at the end of each segment 

c 

IF(ISEG.GT.l) ETAl ( ISEG , IEDGE ) - ETA2 ( ISEG-1 , IEDGE ) 

ETA2 (ISEG, IEDGE) « ETAMX( ISEG, IEDGE) 

C 

c VECTORS FROM ORIGIN TO END POINTS 

C 

DO 80 J — 1 , 3 

RI(J, ISEG, IEDGE) - PT1 (J, ISEG, IEDGE) - COEFE(J, ISEG, IEDGE) 
80 R2(J, ISEG, IEDGE) ■ PT2(J , ISEG, IEDGE) - COEFE(J , ISEG, IEDGE) 
C 

CALL VMAG( Rl(l, ISEG, IEDGE) ,RMAGl) 

C 

CALL VMAG (R2(l, ISEG, IEDGE ),RMAG2) 



C CHOOSE MAPPING FUNCTION 



MAP - MAPSEG( ISEG, IEDGE) 


SEGMAX(J,NMBRSEG( IEDGE) , IEDGE) 
PT2(J, ISEG, IEDGE) 


DO 60 3 — 1 , 3 


ELSE 


PT1(J, ISEG, IEDGE) 
PT2(J, ISEG, IEDGE) 


POINT (J,IPT2( IEDGE) ) 
SEGMAX ( J , ISEG , IEDGE ) 


SEGMAXt J , ( ISEG-1) , IEDGE) 
SEGMAX(J, ISEG, IEDGE) 


END IF 


GO TO (1000,200,300,400,500,600,700) MAP 


C CIRCULAR ARC 

* 

C SWEEP ANGLE 

C 

200 CALL VDOT(Rl(l, ISEG, IEDGE) ,R2(1, ISEG, IEDGE) ,RlDOTR2) 

C 

ARG - R1DOTR2/ ( RMAG1*RMAG2) 


C 

c 

c 

c- 

c 

c- 


IF(ABS(ARG) .GT.1.0) ARG - ARG/ABS( ARG) 
THETA( ISEG, IEDGE) - ACOS(ARG) 

GO TO 1000 
CONICS 


c PROJECTION ONTO CONIC AXIS 

c 

300 IF(COEFE(7, ISEG, IEDGE). GT. 0.0) THEN 
C 


DO 305 1-1,3 


305 

C 


TEMP 

R1(I, ISEG, IEDGE) 
R2( I, ISEG, IEDGE) 

TEMP « RMAG1 


MAP - 3 


- Rl(I, ISEG, IEDGE) 

- R2 (I, ISEG, IEDGE) 

- TEMP 


RMAG1 - RMAG2 
RMAG2 - TEMP 
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C 

C 

c 

c 

c — 
c 

c 

c 

c 

c— 

c 

c 

c- 

c 


c 

c- 

c 


END IF 

CALL VDOT ( Rl(l , ISEG , IEDGE ) , COEFE ( 4 , ISEG , IEDGE ) , Al ) 

CALL VDOT( R2 ( 1 , ISEG, IEDGE) ,COEFE( 4 , ISEG# IEDGE) ,A2) 

—RECIPROCAL OF ECCENTRICITY (RE - 1/E) 

ECT - ABS( ( A1-A2 )/ ( RMAG1-RMAG2) ) 

IF(ECT. GT. 0.9999. AND. ECT. LT. 1.0001) ECT ■ 1.0 
RE (ISEG, IEDGE) « ECT 
—DISTANCE FROM DIRECTRIX TO FOCUS 

RC( ISEG, IEDGE) * RMAG1 * RE ( ISEG, IEDGE) - Al 
— THETA 1 s ANGLE BETWEEN NEGATIVE OF CONIC AXIS AND Rl 

ARG « Al/RMAGl 

I F ( ABS ( ARG ) .GT.1.0) ARG - ARG/ ABS ( ARG ) 

THETAK ISEG, IEDGE) - ACOS(-ARG) 

— THETA2: ANGLE BETWEEN NEGATIVE OF CONIC AXIS AND R2 


C 

C- 

c 

c 

c« 

c 


ARG « A2/RMAG2 

IF(ABS(ARG) .GT.1.0) ARG - ARG/ ABS ( ARG ) 

THETA2 - ACOS(-ARG) 

• — SWEEP ANGLE: THETA2 - THETAl 

THETA ( ISEG, IEDGE) * THETA2 - THETAl ( ISEG, IEDGE) 

LOCAL COORDINATE SYSTEM 

DO 310 J«1 , 3 

UJ(J, ISEG, IEDGE) » Rl (J , ISEG, IEDGE J/RMAG1 
310 UI(J, ISEG, IEDGE) • R2(J,ISEG,IEDGE)/RMAG2 

CALL CROSS (UI( 1, ISEG, IEDGE) ,UJ( 1 , ISEG, IEDGE) ,UK( 1 , ISEG, IEDGE ) ,10) 
INTEGRATION CONSTANTS 

RA( ISEG, IEDGE) - SQRT( ABS (RE (ISEG, IEDGE)** 2 - 1.0)) 

T « TAN (THETAl ( ISEG, IEDGE )*0 . 5) 

— ELLIPSE — — ~ 

IF ( RE (ISEG, IEDGE ) .LT. 1.0001) GO TO 320 

—ARC LENGTH INTEGRAL (CRC EQN. 341) EVALUATED AT THETAl 

ARG - RA ( ISEG, IEDGE )*T/( RE (ISEG, IEDGE) + 1.0) 

ARCH ISEG, IEDGE) - 2.*RC( ISEG, IEDGE )*ATAN( ARG )/RA( ISEG, IEDGE) 

ARC LENGTH INTEGRAL (CRC EQN. 341) EVALUATED AT THETA2 
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T - TAN(THETA2*0.5) 

ARG - RA( ISEG, IEDGE )*T/( RE (ISEG, IEDGE) ♦ 1.0) 

ARC 2 ■ 2.*RC( ISEG, IEDGE)*ATAN(ARG)/RA( ISEG, IEDGE) 

— TOTAL ARC LENGTH 

ARC( ISEG, IEDGE) » ARC 2 - ARC1 ( ISEG , IEDGE ) 

— TANGENT CONSTANT 

XLENGTH( ISEG, IEDGE) • ARC( ISEG, IEDGE)* (RE( ISEG, IEDGE) + 1.0) 
j, /( THETA( ISEG, IEDGE )*RC( ISEG, IEDGE) ) 


GO TO 1000 

— HYPERBOLA ’ 

320 IF ( RE (ISEG, IEDGE ) .GT. 0.9999) GO TO 330 

— ARC LENGTH INTEGRAL EVALUATED AT THETA1 

ARG - (RA( ISEG, IEDGE) *T + RE ( ISEG, IEDGE) + 1.) 
S /( RA( ISEG, IEDGE )*T - RE( ISEG, IEDGE) - 1.) 

ARG - ABS(ARG) 

ARCK ISEG, IEDGE) - RC( ISEG, IEDGE) *LOG(ARG)/RA( ISEG, IEDGE) 

— ARC LENGTH INTEGRAL EVALUATED AT THETA2 

T ■ TAN ( THETA2 *0.5) 

ARG « (RA( ISEG, IEDGE )*T + RE ( ISEG, IEDGE) ♦ 1.) 
t, /(RA( ISEG, IEDGE )*T - RE ( ISEG, IEDGE) - 1.) 

ARG - ABS(ARG) 

ARC 2 » RC ( ISEG, IEDGE )* LOG ( ARG )/RA( ISEG, IEDGE) 


— TOTAL ARC LENGTH 

ARC( ISEG, IEDGE) * ARC 2 - ARCK ISEG, IEDGE) 

TANGENT CONSTANT 

XLENGTH( ISEG, IEDGE) * -4.0*(1.0 + RE ( ISEG, IEDGE) ) 

6 *ARC( ISEG, IEDGE)/(THETA( ISEG, IEDGE)*RC( ISEG, M)) 


GO TO 1000 

— PARABOLA 

— ARC LENGTH INTEGRAL EVALUATED AT THETAl 
330 ARCK ISEG, IEDGE) - RC( ISEG, IEDGE )*T 

ARC LENGTH INTEGRAL EVALUATED AT THETA2 

T » TAN ( THETA2*0 . 5 ) 


— TOTAL ARC LENGTH 

ARC (ISEG, IEDGE) ■ RC( ISEG, IEDGE )*T - ARCl( ISEG, IEDGE) 

C 
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C TANGENT CONSTANT 

C 

XLENGTH( ISEG, IEDGE) ■ 2. 0*ARC( ISEG, IEDGE) 

(, /( THETA ( ISEG# IEDGE ) *RC( ISEG, IEDGE) ) 

GO TO 1000 

EDGE OF REVOLUTION MAP - 4 

400 CONTINUE 

— UK: NORMALIZED VECTOR ALONG AXIS OF .REVOLUTION 
CALL VMAG (COEFE( 4, ISEG, IEDGE ) ,UKM) 

DO 410 1-1,3 

410 UK(I, ISEG, IEDGE) - COEFE(I + 3, ISEG, IEDGE )/UKM 

— UJ: NORMALIZED VECTOR FROM AXIS TOWARD EDGE 

CALL CROSS ( UK ( 1, ISEG, I EDGE ) ,R1(1, ISEG, IEDGE) ,UJ( 1 , ISEG, IEDGE) , 20) 

— UI: NORMALIZED VECTOR FROM AXIS TO FIRST POINT 

CALL CROSS { UJ ( 1 , ISEG , IEDGE ) , UK ( 1 , ISEG , IEDGE ) , UI ( 1 , ISEG , IEDGE ) , 21 ) 

— Rl: VECTOR FROM AXIS TO FIRST POINT 

CALL VDOT( Rl ( 1 , ISEG, IEDGE) ,UK( 1, ISEG, IEDGE) ,RPAl ( ISEG, IEDGE) ) 

CALL VADD (1.0, Rl(l, ISEG, IEDGE ) ,-RPAl( ISEG, IEDGE) ,UK( 1 , ISEG, IEDGE) , 

£, Rl(l, ISEG, IEDGE) ,VECTER) 

CALL VMAG( Rl ( 1 , ISEG, IEDGE ) ,RMl( ISEG, IEDGE) ) 

R2: VECTOR FROM AXIS TO SECOND POINT 

CALL VDOT ( R2 ( 1 , ISEG , IEDGE ) , UK ( 1 , ISEG , IEDGE ) , RPA2 ( ISEG , IEDGE ) ) 

CALL VADD(1.0,R2(1, ISEG, IEDGE) ,-RPA2( ISEG, IEDGE) ,UK( 1 , ISEG, IEDGE) , 
& R2(l, ISEG, IEDGE ),VECTER) 

CALL VMAG( R2 ( 1 , ISEG , IEDGE ) , RM2 ( ISEG , IEDGE ) ) 

SWEEP ANGLE 

CALL VDOT (Rl(l, ISEG, IEDGE) ,R2( 1 , ISEG, IEDGE) ,RlDOTR2) 

ARG - R1DOTR2/ (RM1{ ISEG, IEDGE )*RM2( ISEG, IEDGE ) ) 

I F ( ABS ( ARG ) . GT . 1 • ) ARG - ARG/ABS( ARG) 

THETA( ISEG, IEDGE) - ACOS(ARG) 

GO TO 1000 

USER SUPPLIED SPECIAL EDGE INITIALIZATION MAP - 5 

500 CONTINUE 
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USER SUPPLIED SPECIAL EDGE INITIALIZATION MAP - 6 

600 CONTINUE 

USER SUPPLIED SPECIAL EDGE INITIALIZATION MAP * 7 

700 CONTINUE 

1000 CONTINUE 

RETURN 

END 

•••••mapping**** ******************************************************* 


SUBROUTINE EMAP{ I EDGE, ISEG, RATIO, POINT, TANGENT) 


THIS ROUTINE CALCULATES THE COORDINATES, AND TANGENT 
OF A POINT ON AN EDGE. 


COMMON 

COMMON 

6 

COMMON 

& 

COMMON 

COMMON 

COMMON 


/COEFF/ COEFE ( 8 , 1 0 , 1 2 ) , COEFS ( 8 , 6 ) , NMBRS EG ( 1 2 ) 

/EDGE0/ UI (3, 10, 12), UJ (3, 10, 12), UK (3, 10, 12), 

Rl(3,10,12) ,R2(3,10,12) ,THETA( 10 , 12 ) 

/EDGE3/ ARC( 10,12) ,ARC1 (10,12) ,XLENGTH( 10, 12) , 

RA(10,12),RC(10,12),RE(10,12) ,THETAl (10,12) 
/EDGE8/ RM1{10,12) , RM2(10,12) ,RPA1(10,12) ,RPA2(10,12) 
/MAPING/ MAPS IDE ( 6 ) ,MAPSEG( 10,12) 

/SEGMENT/ PT1(6,10,12) ,PT2(6,10,12) , ETA1(10,12) , ETA2(10,12) 


COMMON /DFN1/ 
COMMON /DFN2/ 
COMMON /DFN3/ 
COMMON /DFN5/ 


DU 1 ( 3 ) , DU2 ( 3 ) , DU3 ( 3 ) , DFNB , DFND, DFNF 

AE , BE , DFNR , ZSTAR , AGL , BETAl , BETA2 , BETA3 , BETA4 

AD,BD,BETA( 37) ,D£LRHO( 37 ) ,RADOC( 37 ) 

CE , CO , ABOT , ATOP 


DIMENS ION XN( 3 ) , XR( 3 ) ,XP( 3 ) , VDUM ( 3 ) 


DIMENSION POINT ( 6 ) , TANGENT ( 3 ) 
DIMENSION VECTOR { 3 ) 


MAP - MAPS EG ( ISEG, IEDGE) 


GO TO (100,200,300,400,500,600,700) MAP 


LINEAR EDGE MAP - 1 


100 DO 110 1-1,3 

POINT(I) - (1. - RATIO) *R1( I, ISEG, IEDGE) + RATIO* R2( I , ISEG, IEDGE) 

110 VECTOR(I) « R2( I, ISEG, IEDGE) - Rl( I, ISEG, IEDGE) 

C 

CALL VMAG( VECTOR, RMAG) 

C 

IF ( RMAG. EQ. 0.0) RMAG « 1.0 
C 

DO 120 1-1,3 

120 TANGENT ( I ) « (R2( I, ISEG, IEDGE) - R1 ( I , ISEG , IEDGE )) /RMAG 
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GO TO 1000 


CIRCULAR ARC MAP » 2 


200 PHI1 - (1.0 - RATIO)*THETA(ISEG, IEDGE) 

PHI2 * RATIO*THETA( ISEG , I EDGE ) 

DO 210 1-1,3 

POINT(I) ■ C0EFE( I , ISEG, I EDGE) + 

1 (SIN(PHI1)*R1(I, ISEG, I EDGE) + SIN ( PHI 2 ) *R2 ( I , ISEG, IEDGE) ) 

2 /SIN(THETA( ISEG, IEDGE)) 

210 TANGENT ( I ) - THETA( ISEG , IEDGE ) * 

1 (COS(PHI2) *R2( I , ISEG, IEDGE) - COS ( PHIl )*Rl ( I , ISEG, IEDGE) ) 

2 /SIM(THETA( ISEG, IEDGE) ) 

GO TO 1000 


CONICS MAP ■ 3 


300 DARC - RATIO* ARC ( ISEG, IEDGE) 

IF(COEFE(7, ISEG, IEDGE) .GT. 0.0) DARC - (1. - RATIO) *ARC( ISEG, IEDGE ) 


— RE • 1/E 


IF(RE( ISEG, IEDGE) .LT. 0.9999) GO TO 310 
IF(RE(ISEG, IEDGE). GT. 1.0001) GO TO 320 

GO TO 330 


— HYPERBOLA 

— CALCULATE ANGLE CORRESPONDING TO ARC LENGTH 

310 ARG1 - (DARC + ARC1( ISEG, IEDGE ) )*RA( ISEG, IEDGE) 
& /RC( ISEG, IEDGE) 

ARG2 - (EXP(ARGl) - 1 . 0 )*( RE ( ISEG , IEDGE ) + 1.0) 
6 /(RA( ISEG, IEDGE)* (EXP (ARG1) + 1.0)) 

ALPHA - 2 . 0*ATAN ( ARG2 ) 

DERIVATIVE 

DEDN - XLENGTH ( ISEG , IEDGE ) *EXP ( ARG1 ) 

& /( (1.0 + ARG2**2)*(1.0 - EXP( ARG1 ) ) **2 ) 

GO TO 340 



C 

C 


—ELLIPSE 

—CALCULATE ANGLE CORRESPONDING TO ARC LENGTH 

320 ARG1 • (DARC + ARC1 ( ISEG, IEDGE )) *RA( ISEG, IEDGE ) *0 . 5 
& /RC( ISEG, IEDGE) 

ARG2 - (RE( ISEG, IEDGE) + 1 . 0 )*TAN(ARG1 )/RA( ISEG, IEDGE ) 

ALPHA > 2 . 0*ATAN ( ARG2 ) 

— DERIVATIVE 

DEDN » XLENGTH( ISEG, IEDGE) *(1.0 + TAN(ARG1)**2)/(1. + ARG2**2) 
IF(COEFE(7, ISEG, IEDGE). GT. 0.0) DEDN - -DEDN 
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GO TO 340 


C PARABOLA — 

C 

c CALCULATE ANGLE CORRESPONDING TO ARC LENGTH 

C 330 ARG - -(DARC + ARC1 ( ISEG, IEDGE) )/RC( ISEG, IEDGE) 

ALPHA - 2 • 0 * ATAN ( ARG ) 

C 

C DERIVATIVE 

C 

DEDN « XLENGTH (ISEG, IEDGE )/( 1*0 + ARG**2) 

C 



C 

C MAGNITUDE OF POSITION VECTOR 

C 340 R « RC( ISEG# IEDGE )/( RE ( ISEG* IEDGE) + COS ( ALPHA) ) 

C 

c RATIO USING ANGLES 

ESP - (ALPHA - THETAl ( ISEG* IEDGE) )/THETA( ISEG* IEDGE) 

C 

PHIl « (1.0 - ESP) *THETA( ISEG, IEDGE) 

PH I 2 - ESP*THETA( ISEG, IEDGE) 

C 

C POSITION AND TANGENT 

C 

DO 350 1-1,3 

POINT ( I ) • COEFE( I, ISEG, IEDGE) + 

1 R* ( SIN ( PHIl )*UJ( I, ISEG, IEDGE) ♦ SIN( PHI2) *UI( I , ISEG, IEDGE) ) 

2 /SIN(THETA(ISEG, IEDGE) ) 


TANGENT ( I ) - (POINT(I) - COEFE( I , ISEG, IEDGE ) ) 

1 * THETA ( ISEG, IEDGE )*R*SIN( ALPHA)/RC( ISEG, IEDGE) 

2 + THETA( ISEG, IEDGE )*R 

3 *( -COS ( PHI 1)*UJ( I, ISEG, IEDGE) + COS( PHI2)*UI ( I , ISEG, IEDGE) ) 

4 /SIN(THETA( ISEG, IEDGE) ) 


C 

C 


350 TANGENT ( I ) 
GO TO 1000 


TANGENT ( I)*DEDN 


MAP - 4 


C EDGE OF REVOLUTION 



c PROJECTION ALONG AXIS AND RADIUS 

400 RP - ( RPA2( ISEG, IEDGE)-RPA1( ISEG, IEDGE) )*RATIO ♦ RPAl ( ISEG, IEDGE) 
RM - ( RM2( ISEG, IEDGE)- RMl ( ISEG, IEDGE) ) ‘RATIO ♦ RM1 ( ISEG, IEDGE ) 

C 

C ANGLE 

C 

GAMMA - THETA( ISEG, IEDGE) ‘RATIO 
C 

C CALCULATE THE POSITION AND TANGENT 

C 

DO 410 1-1,3 
C 

UR - COS ( GAMMA )*UI( I, ISEG, IEDGE) 

& + SIN(GAMMA)*UJ(I, ISEG, IEDGE) 
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C 

C 


POINT( I) 

410 TANGENT ( I ) 
& 


COEFE( I , ISEG, IEDGE) + RP*UK( I, ISEG, I EDGE) + RM*UR 

COS ( GAMMA ) * U J ( I , ISEG , I EDGE ) 

SIN( GAMMA) *UI( I, ISEG, IEDGE) 


C 

GO TO 1000 



C HGM HOLE (EDGES 1,6, 9, 5) 



c ANGULAR LOCATION (B) 

c 

500 BETA I - COEFE( 1,1, IEDGE) 

BETAF - COEFE( 2,1, IEDGE) 


C 


DBTA * BETAF - BETAI 

B - BETAI + RATIO* DBTA 


5 


C 

c CALCULATE HOLE RADIUS (RHO) 

C 

CALL DELRAD ( B , RHO , DRHO , RHOD , RDOC ) 

C 

C CALCULATE AXIAL DISTANCE (CP) 

C 

CALL CAXIS(B,RHO,CP) 

C 

C CALCULATE COORDINATES OF POINT ON HOLE 

C 

CALL HOLE( B , POINT, VDUM) 

C 

C TANGENT 

C 

CB » COS(B) 

SB - SIN(B) 

C 

C RADIAL UNIT VECTOR 

C XR ( 1 ) * CB*DU2( 1 ) + SB*DU3( 1 ) 

XR( 2 ) - CB*DU2 ( 2 ) + SB*DU3(2) 
XR( 3 ) - CB*DU2( 3) + SB*DU3{3) 

C 

CALL DERI V ( B , RHO , CP , DRDB , DCDB ) 


+ TANGENT ( 1 ) ) *DBTA 
♦ TANGBNT( 2 ) ) *D8TA 
+ TANGENT ( 3 ) ) *DBTA 

C 

GO TO 1000 

C- - 

C DUCT EXIT PLANE (EDGES 3,7,8,11) *** 



C ANGULAR LOCATION (B) 

C 

600 BETAI - COEFE( 1,1, IEDGE) 

BETAF * COEFE( 2,1, IEDGE) 

C 


CALL VADD( DCDB , DUl , DRDB , XR, TANGENT , VDUM ) 

C 

C TANGENT 

C 

TANGENT ( 1 ) - ( RHO* (-SB*DU2 ( 1 ) +CB*DU3(1)) 
TANGENT ( 2 ) - (RHO*(-SB*DU2(2) + CB*DU3(2)) 
TANGENT ( 3 ) - ( RHO* (-SB*DU2 ( 3 ) + CB*DU3(3)) 
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DBTA - BETAF - BETAI 

B - BETAI + RATIO*DBTA 

— CALCULATE COORDINATES OF POINTS ON THE EXIT PLANE 
CALL DEXIT(B, POINT, VDUM) 

— TANGENT — — 


CB « COS ( B ) 

SB - SIN(B) 

— DUCT RADIUS (RHOD) 

EX1 « (AD*SB)**2 + (BD*CB)**2 
EX2 » BD**2 - AD* *2 

RHOD » AD*BD/SQRT ( EXl ) 

DRDDB - RHOD*SB*CB*EX2/EXl 
DCDB « (RHOD*CB + SB* DRDDB ) *DFNB/ ( CO - CE) 

CALL VADD ( CB , DU 2 , SB , DU 3 , XR , VDUM ) 

CALL VADD( DCDB , DU1 , DRDDB , XR, XP , VDUM ) 

CALL VADD(-SB,DU2 , CB, DU3 ,XN,VDUM) 

— TANGENT 

TANGE NT ( 1 ) « (XP(1) + RHOD*XN( 1 ) ) *DBTA 
TANGENT ( 2 ) « (XP(2) + RHOD*XN( 2) )*DBTA 
TANGENT ( 3 ) « (XP(3) + RHOD*XN ( 3 ) ) *DBTA 

GO TO 1000 


WELD (EDGES 2,4,10,12) MAP * 7 


ANGULAR LOCATION (BB) 

700 PI2 » 1.57079633 

BB - COEFE (1,1, IEDGE ) 

IF (RATIO.LT. 0.0) BB - POINT(l) 

RATIO - ABS( RATIO) 

— WELD RADIUS (RDOC) & DIFFERENCE BETWEEN HOLE AND DUCT RADIUS ( DRHO ) 

CALL DELRAD(BB,RHO, DRHO, RHOD, RDOC) 

— CALCULATE AXIAL DISTANCE (CP) 

CALL CAXIS(BB,RHO,CP) 

C 

y « RDOC - DRHO 

c 

AAD • 0.5*( (ABOT + ATOP) + (ATOP - ABOT)*COS(BB - PI2)) 

X « SORT ( DRHO* ( 2 . 0*RDOC - DRHO)) 


C-27 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 


noon on oonno non nooononoooooo no 


LMSC-HEC TR F225968 


WELD ANGLE 

TH2 - ATAN ( X/Y ) 

THl - PI2 - TH2 

XX * AAD - X 

TH3 « ATAN ( XX/RDOC ) 

ANGL - RATIO* (TH2 + TH3) + THl 

— POINT 6 TANGENT 

IF( ANGL.GE.PI2) THEN 

— DUCT 


AXIAL DISTANCE 

RADIUS 

COORDINATES 

TANGENT 


WELD' 


PSI - (ANGL - PI2)/TH3 
CS - CP + X + PSI*XX 
RS ■ RHOD 

CALL DUCT (BB,RS,CS, POINT) 

TANGENT ( 1 ) - DUl(l) 

TANGENT ( 2 ) » DUl(2) 

TANGENT ( 3 ) » DU1(3) 

RETURN 

ELSE 


XS - RDOC * COS ( ANGL ) 

YS - RDOC* SIN (ANGL) 

— AXIAL DISTANCE 

CS - CP + X - XS 

— RADIUS 

RS - RHOD + RDOC - YS 

— COORDINATES 

CALL DUCT (BB,RS,CS, POINT) 

— TANGENT 

DXDN - TAN (ANGL) 

CB - COS(BB) 

SB - SIN(BB) 

CALL VADD ( CB , DU 2 , SB , DU 3 , VDUM , TANG ENT ) 
CALL VADD ( DXDN , DU1 , - 1 . 0 , VDUM , C , TANGENT ) 
RETURN 
END IF 


1000 CONTINUE 
C 

RETURN 

C 
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1100 DO 1110 1-4,6 

1110 POINT ( I ) - PT1 ( I , ISEG, IEDGE) 

RETURN 

END 


♦•"MAPPING 


SUBROUTINE SURFACE ( INIT, ISIDE) 


THIS ROUTINE DETERMINES THE COORDINATES AND NORMAL OF A 
POINT ON A SURFACE. 


COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


/COEFF/ COEFE( 8 , 10 , 12 ) ,COEFS(8,6) , NMBRSEG(12) 

/COORD/ UA( 3 , 6 ) , UE ( 3 , 6 ) , UN ( 3 , 6 ) 

/COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
/HEADER/ ITITLE( 20) , LINE 
/INITA/ I Z INDEX, MAPTEN, INCHES 

/INITB/ IBOX (4,6), IRAY (4,3), NRAY , IPT1 ( 1 2 ) , IPT2 ( 1 2 ) 

/INITC/ PI ,RADDEG 

/INPUTA/ EDGE( 3,12) ,POINT( 3,8) ,SIDE( 3,6) 

/MARCHS/ MARCH,INDEX(3) 

/MAPED/ KSEG (12), UAXIS ( 3,6), IEDGE1 ( 6 ) , IEDGE 2 ( 6 ) , GAMMA 

/MAPING/ MAPSIDE(6) ,MAPSEG( 10,12) 

/MAXIMUM/ ETAMAX(10,12) ,SEGMAX( 3, 10, 12) 

/PARTIAL/ DEDN (3,12), DSDN (3,2), SNORMAL (3,6) 

/SPACING/ ISTRTCH ( 3 ) , STRETCH ( 3 ) , ETAS (3,200), ETA( 3 ) , DETA( 3 ) 
/UNITS/ NU5,NU6,NU20 

/ZONING/ IZONE , ISECT , NZ INDEX , NMBRNDS ( 3 ) 


COMMON /DFNl/ 
COMMON /DFN2/ 
COMMON /DFN3/ 
COMMON /DFN4/ 
COMMON /DFN5/ 


DUl ( 3) ,DU2( 3 ) , DU3( 3 ) , DFNB, DFND,DFNF 

AE , BE , DFNR , ZSTAR , AGL , BETAl , BETA2 , BETA3 , BETA 4 

AD,BD,BETA( 37) , DELRHO(37) ,RADOC( 37) 

IHOLE,PL( 6, 6, 4), PA (2, 6, 4) 

CE , CO , ABOT , ATOP 


DIMENSION RAM ( 3 ) , SAC ( 3 ) 

DIMENSION Rl(6) ,R2(6),R3<6) ,R4(6) 
DIMENSION E1(6),E2(6),E3(6),E4(6) 


DIMENSION F( 8) 

DIMENSION NMBRPT( 4 ,6) ,SPOINT( 3 ,4 , 6 ) , SEDGE ( 3 , 4 ) , SWEEP ( 3 ) 
DIMENSION UNORMAL( 3,6) ,URADIAL( 3,6) 

DIMENSION ZERO( 3 ) , VECTOR( 3 ) 

DIMENSION IETA1(6) ,IETA2(6) ,SIGNS(6) 

DIMENSION UNITRAD( 3 ) 


DATA IETA1 /I , 1 , 1 , 1 , 2, 2/ 

DATA IETA2 /2, 3, 2, 3, 3., 3/ 

DATA SIGNS /I . 0 ,-l . 0 ,-l . 0 , 1 . 0 , 1 . 0 ,-l . 0/ 
DATA ZERO /0.,0.,0./ 


NODNUM - NODESAV + NODETOT + 1 


— ETA COEFFICIENTS FOR BI-LINEAR INTERPOLATION 


F(l) - 1.0 - ETA( IETA1( ISIDE) ) 
F( 2) - ETA( IETA1( ISIDE) ) 


C-29 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 



LMSC-HEC TR F225968 


F( 3 ) - 1.0 - ETA( IETA2( ISIDE) ) 

F ( 4 ) » ETA(IETA2( ISIDE)) 

F(5) - F( 1 ) *F{ 3 ) 

F( 6 ) ■ F ( 1 ) *F( 4 ) 

F(7) - F ( 2 ) *F( 3) 

F( 8 ) * F(2)*F<4) 

C 

C DETERMINE THE EDGES OF THE SURFACE 

C 

LINE1 « IBOX(l, ISIDE) 

LINE2 « IBOX( 2 , ISIDE ) 

LINE3 - IBOX( 3 , ISIDE) 

LINE4 - IBOX( 4 , ISIDE) 

C 

C DETERMINE THE CORNER POINTS OF THE SURFACE 

C 

LPT1 - IPT1( LINED 
LPT2 - IPTKLINE3) 

LPT 3 ■ IPT2(LINE3) 

LPT 4 - IPT2( LINED 

C 

MAP * MAPSIDE( ISIDE) 


C 

GO TO (100,100,300,400,500,600,700) MAP 



C FLAT OR CYLINDRICAL SURFACE MAP » 1 OR 2 


C 

C CALCULATE POSITION 

C 

100 DO 110 J«l,3 
C 

SIDE(J, ISIDE) 

1 

2 

3 


« F(1)*EDGE(J, LINED 
+ F ( 2 ) *EDGE ( J , LINE3 ) 

♦ F ( 3 ) *EDGE( J , LINE2 ) 

♦ F( 4 ) *EDGE( J , LINE4 ) 


- F( 5 ) *POINT{ J, LPT1 ) 

- F{ 6 ) *POINT( J, LPT4 ) 

- F { 7 ) * POINT { J , LPT2 ) 

- F( 8 ) *POINT( J, LPT3 ) 


C 


C 


C 

c 


1 

2 

3 


DSDN ( J , 1) 


EDGE( J , LINE3 ) 
• EDGE(J, LINED 
+ F ( 3 ) * DEDN ( J , LINE 2 ) 
+ F( 4 ) *DEDN( J , LINE4 ) 


+ F(3)*POINT(J,LPTl) 
+ F(4)*POINT(J,LPT4) 

- F(3)*POINT(J,LPT2) 

- F( 4 ) * POINT (J , LPT3 ) 


no 


1 

2 
3 


DSDN ( J , 2 ) 


EDGE(J,LINE4) 
- EDGE(J,LINE2) 
+ F(1)*DEDN(J, LINED 
+ F ( 2 ) * DEDN ( J , LINE3 ) 


+ F( 1) * POINT (J,LPT1) 
+ F ( 2 ) *POINT ( J , LPT2 ) 

- F( 1) * POINT (J , LPT4) 

- F( 2 ) *POINT( J , LPT3 ) 


CALL CROSS ( DSDN (1,1), DSDN (1,2) ,SNORMAL( 1 , ISIDE) , 41) 


GO TO 1000 




C SPECIAL SURFACE MAP » 3 



300 GO TO 1000 


C EDGE OF REVOLUTION MAP - 4 



400 IF(INIT.EQ.l) THEN 
C 

C ETA DIRECTION OF REVOLUTION 

C 
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IONETWO ■ 1 

I F ( COEFS ( 7 , IS IDE) . EQ. IETA2 ( ISIDE ) ) IONETWO ■ 2 
EDGE 1 

IEDGEK ISIDE) - IBOX( IONETWO, ISIDE) 

EDGE 2 

IEDGE2( ISIDE) ■ IBOX( IONETWO + 2, ISIDE) 

UA: NORMALIZED VECTOR ALONG AXIS OF REVOLUTION 

CALL VADD(1. 0,COEFS(4, ISIDE) ,1.0, ZERO, VECTOR, UA(1, ISIDE) ) 

UN: NORMALIZED VECTOR FROM AXIS TOWARD SURFACE 

CALL VADD(1.0,EDGE(1, IEDGEK ISIDE) ) ,-1.0,COEFS( 1, ISIDE) , El, VECTOR) 
CALL CROSS(UA(l, ISIDE) ,E1 ,UN( 1 , ISIDE ) , 42 ) 

UE: NORMALIZED VECTOR FROM AXIS TO FIRST EDGE 

CALL CROSS ( UN ( 1, ISIDE ) ,UA(1, ISIDE ) ,UE(1, ISIDE) ,43) 

END IF 

Rl: VECTOR FROM AXIS TO FIRST EDGE 

CALL VADD( 1.0, EDGE ( 1 , IEDGE1 ( IS IDE ) ) , -1 . 0 , COEFS ( 1 , ISIDE ) , El , VECTOR ) 
CALL VDOTtEl, UA(1, ISIDE) ,EPA1) 

CALL VADD( 1.0, El, -EPA1,UA(1, ISIDE ) ,R1, VECTOR) 

CALL VMAG(R1,RM1) 

DERIVATIVE 

CALL VDOT(DEDN(l, IEDGEK ISIDE) ) ,UA(1, ISIDE) ,DA1) 

CALL VDOT(DBDN( 1, IEDGEK ISIDE) ) , VECTOR, DR1) 

CALL VMAG(DEDN( 1, IEDGEK ISIDE) ) ,DM1 ) 

R2: VECTOR FROM AXIS TO SECOND EDGE 

CALL VADD( 1.0, EDGE ( 1 , I EDGE 2 ( ISIDE)), -1.0, COEFS ( 1 , ISIDE ) , E2 , VECTOR) 
CALL VDOT(E2,UA(l, ISIDE), EPA2) 

CALL VADD(1.0,E2,-EPA2,UA(1, ISIDE) ,R2, VECTOR) 

CALL VMAG(R2,RM2) 

DERIVATIVE 

CALL VDOT ( DEDN ( 1 , IEDGE2 ( ISIDE) ) ,UA( 1 , ISIDE) ,DA2) 

CALL VDOT ( DEDN ( 1 , IEDGE2( ISIDE) ) , VECTOR, DR2) 
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C 

CALL VMAG(DEDN( 1 , IEDGE2( ISIDE) ) ,DM2) 

— SWEEP ANGLE 

CALL VDOT ( R1 , R2 , R1D0TR2 ) 

ARG - R1D0TR2/ (RM1*RM2) 
I F ( ABS ( ARG ) . GT .1.0) ARG - ARG/ABS ( ARG ) 

THETA « ACOS ( ARG ) 

— CALCULATE RADIUS 

N - COEFS ( 7 , ISIDE ) 

— POSITION 

EP - (EPA2 - EPA1 )*ETA(N) + EPA1 
RM - ( RM2 - RM1 )*ETA(N) + RM1 

— DERIVATIVE 

DA - (DA2 - DA1)*ETA(N) + DAI 
DR - (DR2 - DRl ) *ETA(N) + DR1 
DM « (DM2 - DM1 )*ETA(N) + DM1 


— ANGLE 

GAMMA « THETA* ETA(N) 

— CALCULATE THE POSITION AND SURFACE NORMAL' 
DO 410 1-1,3 


UR - 

SIDE(I, ISIDE) - 
DSDN (1,1) - 
410 DSDN (1,2) - 


COS(GAMMA)*UE(I, ISIDE) + SIN ( GAMMA) *UN< I , ISIDE) 
COEFS (I, ISIDE) + £P*UA( I , ISIDE ) RM*UR 

(DR*UR + DA*UA( I , ISIDE) )/DM 

COS ( GAMMA) *UN( I, ISIDE) - SIN ( GAMMA) *UE( I, ISIDE) 


-ORIENT SURFACE NORMAL DEPENDING ON DIRECTION OF EDGE REVOLUTION 
IF ( COEFS (7, ISIDE ).E0.IETA2( ISIDE) ) THEN 
CALL CROSS ( DSDN (1,1), DSDN (1,2), SNORMAL ( 1 , IS IDE ) , 44 ) 

ELSE 

CALL CROSS ( DSDN ( 1 , 2 ) , DSDN (1,1), SNORMAL ( 1 , ISIDE ) , 45 ) 

END IF 


GO TO 1000 

BOWL SURFACE WITH HOLE MAP » 5 

500 I HOLE - 0 
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IET1 » 1 

IF(ETA( 1 ) .GE .ETAMAX( 1 » 3 ) ) IET1 » 2 
IF(ETA(1) .GE.ETAMAX(2,3) ) IET1 - 3 
IF(ETAU).GE.ETAMAX(3,3)) IET1 - 4 
I F ( ETA ( 1 ) . GE . ETAMAX (4,3)) IET1 ■ 5 

STR11 - 2.0 
STR12 -6.0 
STR13 - 2.0 

EPS1 » ETA( 1 )/ETAMAX( 1 » 3 ) 

IF ( IET1 . EQ. 2 .OR. IET1.EQ.3) THEN 

RATIOl - (ETA(l) - ETAMAX ( 1 r 3 ) ) / (ETAMAX(3,3) - ETAMAX(1,3)) 
XI - RATI01*STR11/2.0 

ETAX1 * ( EXP ( XI ) - EXP ( -XI) )/(EXP(Xl ) + EXP(-Xl)) 

X2 - STR11/2 . 0 

ETAX2 » (EXP(X2) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 

EPX1 « ETAX1/ETAX2 


END IF 


IF ( IET1.EQ.4) THEN 

RATIOl - (ETA(l) - ETAMAX( 3 » 3) ) / (ETAMAX(4,3) - ETAMAX(3,3)) 
XI - (STR12/2.0)/2.0 

ETAMID « ( EXP (XI ) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

X2 - (RATIOl - 0.5)*STR12/2.0 
ETAMX « ( EXP ( X2 ) - EXP(-X2) )/(EXP(X2) ♦ EXP(-X2)) 

EPS1 - (ETAMX + ETAMID)/ ( 2. 0*ETAMID) 

END IF 


IF(IETl.EQ.S) THEN 

RATIOl « (ETA(l) - ETAMAX ( 4,3)) / (1.0 - ETAMAX(4,3)) 
XI ■ (STR13/2.0)/2.0 

ETAMID - (EXP(Xl) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

X2 * (RATIOl - 0. 5) # STR13/2.0 
ETAMX « ( EXP ( X2 ) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 

EPS1 - (ETAMX + ETAMID) /( 2. 0*ETAMID) 

END IF 


IET3 ■ 1 

IF(£TA(3) .GE.ETAMAX(1*8) ) IET3 ■ 2 
IF( ETA( 3 ) .GT.ETAMAX( 2 , 8 ) ) IET3 ■ 3 

STR31 -2.0 
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STR32 - 4.0 
STR33 « 2.0 

IF(IET3.EQ.l) THEN 

RATI03 - ETA(3)/ETAMAX(1,8) 

XI - RAT 103* STR3 1/2.0 

ETAX1 » ( EXP (XI ) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

X2 - STR31/2.0 

ETAX2 - (EXP(X2) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 

EPS3 - ETAX1/ETAX2 
END IF 

IF(IET3.EQ.2) THEN 

RATI03 - ( ETA( 3 ) - ETAMAX(1,8) )/(ETAMAX( 2,8) - ETAMAX(1,8)) 
XI » (STR32/2.0)/2.0 

ETAMID « (EXP (XI) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

X2 - ( RATI03 - 0.5)*STR32/2.0 
ETAMX - ( EXP ( X2 ) - EXP(-X2))/(EXP(X2) + EXP(-X2)) 

EPS3 - (ETAMX + ETAMID)/( 2.0*ETAMID) 

END IP 

I F ( IET3.EQ.3) THEN 

RATI03 - 1.0 - (ETA(3) - ETAMAX( 2,8 ) )/( 1 .0 - ETAMAX(2,8)) 
XI » RATIO3*STR33/2.0 

ETAX1 - ( EXP (XI ) - EXP(-X1))/(EXP(X1) ♦ EXP(-Xl)) 

X2 - STR33/2.0 

ETAX2 - ( EXP ( X2 ) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 

EPS 3 - 1.0 - ETAX1/ETAX2 
END IF 

-BOWL AXIS 

UAXIS ( 1 , IS IDE) « 1.0 
UAXIS ( 2 , ISIDE ) « 0.0 
UAXIS ( 3 , ISIDE) ■ 0.0 

-ARC CENTER ON AXIS 

SAC ( 2 ) ■ 0.0 
SACO) ■ 0.0 

DO 510 J«l,6 

-SUB-SURFACE CORNER POINTS 


c-3a 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 



ono non non non non non 9°*? ° ono 


LMSC-HEC TR F225968 



R1(J) » PL( J , IET1 , IET3 ) 

R2( J) - PL(J,IET1+1,IET3) 

R3(J) - PL( J » IET1+1 , IET3+1 ) 

R4(J) - PL(J,IET1 , IET3+1 ) 

— SUB-SURFACE EDGE POINTS 

El( J) » R1(J) + EPS1* (R2( J) - R1(J) ) 

E2(J) - R2(J) + EPS3*(R3(J) - R2(J)) 

E3(J) » R4 ( J ) ♦ EPS1*(R3(J) - R4(J)) 

510 E4(J) « R1(J) + EPS3MR4 ( J) - R1(J)) 

GO TO (520,540,560) IET3 


IET3 - 1 


520 GO TO (522,524,524,528,530) IET1 
— IET1 ■ 1 

522 DO 523 J«l,6 

E1(J) » EDGE( J , 3 ) 

523 E4 ( J ) - EDGE( J ,8) 

SAC ( 1 ) - SEGMAX( 1 , 1 , 3 ) 

CALL SPEDGE(R2,R3,EPS3,SAC,E2) 

GO TO 580 

— IET1 ■ 2 


24 R2(I) - PL( 1 ,4,IET3) 

R3(l) - PL( 1,4, IET3+1 ) 

E 2 ( 1 ) • R2( 1) + EPS3*(R3(1) - R2(l)) 

EXM « SEGMAX( 1,1,3) + EPX1*(E2(1) - SEGMAX( 1 , 1 , 3 ) ) 

IF ( EXM • GT. SEGMAXf 1,2,3)) GO TO 526 

— INTERIOR 


EPS1 - (EXM -SEGMAX( 1,1,3) )/(SEGMAX( 1,2,3) -SEGMAX( 1,1,3)) 


SUB-SURFACE EDGE POINTS 


C 

c 



c 


525 

C 


IET1 » 2 
DO 525 J»l,6 

R1(J) - PL( J , IET1 , IET3) 

R2(J) - PL( J, IETI+1 , IET3) 

R3(J) « PL( J,IET1+1, IET3+1 ) 

R4( J) » PL( J , IET1 , IET3+1 ) 

El ( J) • EDGE(J,3) 

E3 ( J) ■ R4(J) + EPS1*(R3(J) - R4(J)) 
SAC ( 1 ) ■ SEGMAX( 1 , 2, 3) 
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CALL SPEDGE (R2,R3,EPS3, SAC, E2) 

SAC ( 1 ) - SEGMAXt 1,1,3) 

CALL SPEDGE (R1,R4, EPS 3, SAC, E4) 

GO TO 580 

— IET1 » 3 — 

526 EPS1 - (EXM - SEGMAX( 1 , 2, 3) )/( E2( 1 ) - SEGMAX( 1 , 2 , 3 ) ) 
— SUB-SURFACE EDGE POINTS 

IET1 » 3 
DO 527 J«l,6 

R1(J) ■ PL( J , IETl , IET3) 

R2(J) - PL(J,IET1+1,IET3) 

R3( J) ■ PL( J , IET1+1 , IET3+1 ) 

R4(J) - PL( J , IETl , IET3+1 ) 

E1(J) ■ EDGE(J,3) 

527 E3(J) ■ R4(J) + EPS1*(R3(J) - R4(J)) 

SAC ( I ) - SEGMAX( 1,3,3) 

CALL SPEDGE ( R2 , R3 , EPS 3 , SAC > E 2 ) 

SAC ( 1 ) * SEGMAXt 1,2,3) 

CALL SPEDGE (R1,R4,EPS3, SAC, E4) 

GO TO 580 

— IETl » 4 

528 DO 529 J«l,6 

529 El ( J ) ■ EDGE ( J , 3 ) 

SAC ( 1 ) * SEGMAXt 1 #4 , 3) 

CALL SPEDGE (R2,R3,EPS3, SAC, E2) 

SAC ( 1 ) • SEGMAXt 1 , 2 , 3 ) 

CALL SPEDGE (R1,R4, EPS 3, SAC, E4) 

B « BETA1 + EPS1* (BETA2 - BETA1) 

CALL HOLE (B,E 3, VECTOR) 

GO TO 580 

IETl » 5 

530 DO 531 J«l,6 
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C 

C 

c 

c 


E 1 ( J ) - EDGE(J,3) 

531 E2(J) - EDGE(J,7) 

SAC ( 1 ) - SEGMAX(1,4,3) 

CALL SPEDGE ( R1 , R4 , EPS 3 /SAC , E4 ) 

GO TO 580 


IET3 ■ 2 


540 GO TO (542,544,544,548,550) IET1 

— IET1 - 1 

542 DO 543 J-1,6 

543 E4(J) - EDGE( J ,8) 

SAC ( 1 ) - SEGMAXt 1,1,3) 

CALL SPEDGE (R2,R3,EPS3, SAC, E2) 

GO TO 580 

--IET1 - 2— 

544 B • BETA1 + EPS3*(BETA4 - BETA1) 

CALL HOLE (B,E 2, VECTOR) 

EXM - SEGMAX( 1,1,3) + EPX1*(E2(1) - SEGMAX( 1,1,3)) 
IF ( EXM.GT. SEGMAXt 1,2,3) ) GO TO 546 
— INTERIOR 


EPS1 - (EXM -SEGMAX( 1,1,3) )/( SEGMAX( 1,2,3) —SEGMAXt 1,1,3)) 


SUB-SURFACE EDGE POINTS 


IET1 ■ 2 


C 

C 


c 


c 

c 

c 


DO 545 J»l,6 

R1(J) « PL( J , IET1 , IET3) 
R2(J) ■ PL( J, IET1+1 , IET3 ) 
R3(J) ■ PL( J , IET1+1 , IET3+1 ) 
R4(J) - PL( J, IET1 , IET3+1 ) 

E1(J) ■ R1(J) + EPS1*(R2(J) 
545 E3(J) ■ R4( J) + EPS1*(R3(J) 

SAC ( 1 ) » SEGMAX( 1,2,3) 

CALL SPEDGE (R2,R3,EPS3, SAC, E2) 


C 


SAC ( 1 ) • SEGMAX(1,1,3) 


R1(J) ) 
R4(J) ) 
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CALL SPEDGE ( R1 , R4 , EPS3 , SAC , E4 ) 

GO TO 580 

— IET1 « 3 

546 EPS1 « (EXM - SEGMAX( 1,2,3) )/( E2( 1 ) - SEGMAX( 1 , 2 , 3 ) ) 
-SUB-SURFACE EDGE POINTS 

IET1 - 3 
DO 547 J-1,6 

RKJ) - PL( J , IET1 , IET3) 

R2(J) - PL( J , IET1+1 , IET3 ) 

R3(J) ■ PL( J , IET1+1 , IET3+1 ) 

R4(J) - PL( J , IET1 , IET3+1 ) 

E1(J) « RKJ) + EPS1*(R2(J) - Rl(J)) 

547 E3(J) - R4(J) + EPS1*(R3(J) - R4(J)) 

SAC ( 1 ) » SEGMAX(1,2,3) 

CALL SPEDGE(R1,R4,EPS3,SAC,E4) 

B - BETA1 + EPS3*(BETA4 - BETA1 ) 

CALL HOLE(B,E2, VECTOR) 

GO TO 580 

— IET1 » 4 

548 IHOLE * 1 
B ■ BETA1 + EPS1*(BETA2 - BETA1) 

CALL HOLE (B, El, VECTOR) 

B » BETA2 + EPS3*(BETA3 - BETA2 - 2.*PI) 

CALL HOLE(B,E2, VECTOR) 

B • BETA 4 * EPS1* (BETA3 - BETA4 ) 

CALL HOLE (B,E 3, VECTOR) 

B « BETA1 ♦ EPS3* (BETA4 - BETA1 ) 

CALL HOLE(B,E4, VECTOR) 

GO TO 580 

IET1 - 5 

550 DO 551 J-1,6 
C 

551 E2(J) - EDGE( J , 7 ) 

C 
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B - BETA 2 + EPS3*(BETA3 - BETA 2 - 2.*PI) 

C 

CALL H0LE(B,E4, VECTOR) 

C 

GO TO 580 


IET3 - 3 


560 GO TO (562,564,564,568,570) IET1 
— IET1 - 1 

562 DO 563 J-1,6 

E3(J) - EDGE( J , 11 ) 

563 E4(J) - EDGE( J , 8 ) 

SAC ( 1 ) « SEGMAX( 1,1,3) 

CALL SPEDGE ( R2 , R3 , EPS3 , SAC , E2 ) 

GO TO 580 


— IET1 « 2 


564 R2(l) - PL( 1 ,4 , IET3 ) 

R3(l) « PL( 1,4,IET3+1) 

E2( 1 ) - R2( 1) + EPS3* ( R3 ( 1 ) - R2(l)) 

EXM - SEGMAX( 1,1,3) + EPX1*(E2(1) - SEGMAX( 1 , 1 , 3 ) ) 

IF ( EXM .GT.SEGMAX( 1,2,3)) GO TO 566 


INTERIOR 


C 

C 

C 


c 


565 

C 

C 

c 

c 

c 


EPS1 - (EXM -SEGMAX( 1,1,3) )/( SEGMAX( 1,2,3) 
IET1 « 2 
DO 565 J-1,6 

R1(J) - PL( J , IET1 , IET3) 

R2(J) - PL( J, IET1+1 , IET3 ) 

R3 ( J ) - PL( J , IET1+1 , IET3+1 ) 

R4(J) - PL( J , IET1 , IET3+1 ) 

El ( J ) - R1(J) + EPS1* ( R2 ( J) - R1(J) ) 

E3(J) - EDGE( J , 11 ) 

SAC ( 1 ) - SEGMAX( 1,2,3) 

CALL SPEDGE (R2,R3,EPS3, SAC, E2) 

SAC ( 1 ) « SEGMAX( 1 , 1 , 3 ) 

CALL SPEDGE (R1,R4,EPS3, SAC, E4) 


C 


GO TO 580 


SEGMAX( 1,1,3)) 
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C IET1 * 3 


C 

566 

C 


EPS1 - (EXM - SEGNAX( 1,2,3) )/( E2( 1 ) 
IET1 - 3 


C 

DO 567 J»l,6 

C 

R1 ( J ) - PL( J, IET1 ,IET3) 

R2(J) - PL(J,IET1+1,IET3) 

R3( J) » PL( J , IET1+1 f IET3+1 ) 

R4(J) - PL{ J , IET1 , IET3+1 ) 

C 

E1(J) » Rl( J) + EPS1*(R2(J) - R1(J) ) 
E2( J) « R2( J) + EPS3* ( R3 ( J ) - R2(J)) 
567 E3(J) - EDGE( J , 11 ) 

C 

SAC ( 1 ) « SEGMAX( 1,2, 3) 

C 

CALL SPEDGE ( R1 , R4 , EPS 3 , SAC , E4 ) 

C 

GO TO 580 


SEGMAX(1,2,3) ) 


C 

C IET1 » 4 

C 

568 DO 569 J-1,6 
C 

569 E3(J) - R4( J) + EPS1*(R3(J) - R4(J)) 

C 

B - BETA4 + EPS1*(BETA3 - BETA4) 

C 

CALL HOLE (B, El, VECTOR) 

C 

GO TO 580 
C 

C IET1 - 5 

C 

570 DO 571 J«l,6 
C 

E3( J) - EDGE( J , 11 ) 

571 E2(J) « EDGE( J,7) 

C 



C INTERPOLATION 



C POSITION 

c 

580 DO 582 J«l,3 

C 582 RAM(J) - (1. - EPS1)*(E4(J) - (1. - EPS3)*R1(J) - EPS3*R4(J)) 

1 + EPS1*(E2( J) - (1. - EPS3)*R2{J) - EPS3 # R3(J)) ♦ E1(J) 

2 + EPS3*(E3(J) - El ( J ) ) 

C 



C RADIUS AND TANGENT 

C 

C AXIAL DISTANCE 

C 

RAX - RAM(l) 

C 
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C REGION 1 

C 

IF( RAX.GT .SEGMAX( 1,1,3) ) GO TO 586 


C 

C RADIUS 


C 


PSI - 


C 


c 


RADI » 
RAD2 » 

RAD » 


C 

C TANGENT 

C 

DSDN (1,1) • 
DSDN (2,1) » 
DSDN( 3 r I) ■ 
C 

GO TO 590 
C 

C REGION 2 

C 


RAX/SEGMAX (1,1,3) 

SQRT( POINT( 2»4)**2 + POINT( 3 , 4 )**2) 
SQRT( SEGMAX( 2, 1 , 3 ) **2 + SEGMAX(3,1,3)**2) 

RADI + PSI* ( RAD2 - RADI) 


SEGMAX(1,1,3) - POINT ( 1,4) 
RAD2 - RADI 
0 . 


586 IF(RAX.GT.SEGMAX(1,2, 3) ) GO TO 588 


C 


C RADIUS 


C 


PSI 


C 


c 


RADI 

RAD2 

RAD 


C 

C TANGENT 

C 

DSDN (1,1) 
DSDN (2*1) 
DSDN (3,1) 
C 

GO TO 590 
C 

C REGION 3 

C 


(RAX - SEGMAX( 1,1,3) )/( SEGMAXt 1,2,3) - SEGMAX( 1,1,3)) 

SQRT( SEGMAX( 2 , 1 , 3 ) **2 + SEGMAX(3,1,3)**2) 

SQRT(SEGMAX( 2, 2, 3 )**2 ♦ SEGMAX(3,2,3)**2) 

RADI + PSI* (RAD 2 - RADI) 


SEGNAXt 1,2,3) - SEGMAX( 1,1,3) 
RAD2 - RADI 
0 . 


C ELLIPTIC REGION OF BOWL 


C 

O 

c 

c 

o 

c 

c 

c 


c 


— RADIUS 
588 RAD 

— TANGENT 

DRDX 

TH 

DSDN (1,1) 
DSDN (2,1) 
DSDN( 3, 1 ) 


DFNR + BE*SQRT( 1 .0 - ((RAX - DFND)/AE)**2) 

-( (BE/AE)**2)*(RAX - DFND)/(RAD - DFNR) 

ATAN(DRDX) 

COS(TH) 

SIN(TH) 

0.0 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


INSIDE HOLE 

IF(IET3.EQ.2 .AND. IET1.EQ.4) THEN 
1-0 

THETA ■ AT AN ( RAM( 2) / RAM( 3 ) ) 

DO 587 J*1 , 36 

CALL HOLE(BETA(J) , El, VECTOR) 

THETA1 - ATAN( El(2) / El (3) ) 

CALL HOLE(BETA(J+l) ,E2, VECTOR) 

THETA2 ■ ATAN( E2(2) / E2(3) ) 

MATCH ■ 0 

I F ( THETA . GT . THETA1 .AND. THETA. LT.THETA2) MATCH ■ 1 
IF ( THETA. GT.THETA2 .AND. THETA. LT.THETA1 ) MATCH ■ 2 

IF ( MATCH. GT.O) THEN 

I - I + 1 

I F ( MATCH . EQ . 1 ) THEN 

RATIO - (THETA - THETA1) / ( THETA2 - 
B ■ RATIO* (BETA(J+1) - BETA(J)) 

ELSE 

RATIO » (THETA - THETA2) / (THETA1 - 
B ■ RATIO* (BETA(J) - BETA( J+l ) ) 
END IF 


C 

C 

C 

c 

c 

c 

c 

c 

c 


IF(I.EQ.l) THEN 

CALL HOLE (B,E 3, VECTOR) 

RADI - SQRT( E3(2)**2 + E3(3)**2) 
GO TO 587 
ELSE 

CALL HOLE (B,E4, VECTOR) 

RAD2 « SQRT( E4(2)**2 + E4(3)**2) 
GO TO 589 
END IF 

END IF 


587 CONTINUE 
589 CONTINUE 
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C 

c 

c 

c 

c 

c 

c 

c 


IF(E3(1).LT.E4(1)) THEN 

RATIO - (RAM(l) - E3( 1) ) / (E4(l) - E3(l)) 
RAD - RATIO* (RAD2 - RADI) + RADI 
ELSE 

RATIO - (RAM(l) - E4( 1) ) / (E3(l) - E4(l)) 
RAD ■ RATIO* (RADI - RAD2) + RAD2 
END IF 


DSDN( 1 , 1 ) - 1.0 
DSDN( 2,1) » 0.0 
DSDN (3,1) - 0.0 


C BELOW HOLE 

C 

IF( IET3 . EQ. 3 .AND. RAX. LE .SEGMAX( 1 , 4 , 11 ) ) THEN 

C 

RAD ■ SQRT(RAM( 2)**2 + RAM(3)**2) 

C 

DSDN (1,1) - 1.0 
DSDN(2,1) - 0.0 
DSDN (3,1) - 0.0 


END IF 


C- 

c 

c- 

c- 

c 

c 

c- 

c 

c 

c- 

c 

c 

c- 

c 


OUTPUT POSITION AND NORMAL 
— AXIAL COMPONENT OF THE TANGENT 
590 CALL VDOT ( DSDN (1,1)* UAXIS ( 1 , ISIDE ) , DA) 
— NORMAL COMPONENT OF THE TANGENT 
DN - DSDN (2,1) 

— INTERPOLATED RADIUS 

RADX - SORT ( RAM ( 2 ) * * 2 + RAM(3)**2) 
ANGULAR LOCATION 

ANG - AS IN ( RAM ( 3) /RADX) 

C 

CANG - COS (ANG) 

SANG - SIN(ANG) 

C 

C TANGENT 1 (AXIAL DIRECTION) 

C 

DSDN(1,1) » DA 
DSDN (2,1) ■ DN*CANG 
DSDN (3,1) ■ DN*SANG 
C 

C TANGENT 2 (CIRCUMFERENTIAL DIRECTION) 
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C 

DSDN ( 1,2) - 0.0 

DSDN (2,2) - -SANG 
DSDN(3,2) - GANG 
C 

C NORMAL 

C 

CALL CROSS(DSDN(l, 1) ,DSDN(1, 2) , SNORMAL( 1 , ISIDE ) , 45) 


C 

C POSITION 

C 


c 

c 


IF( EPS3. EQ. 1 .0 ) 


C 


SIDE(1, ISIDE) 
SIDE(2, ISIDE) 
SIDE( 3 t ISIDE) 

THEN 

SIDE(1, ISIDE) 
SIDE(2, ISIDE) 
SIDE ( 3 f ISIDE) 

END IF 


RAX 

RAD* GANG 
RAD* SANG 


E3(l) 
E3( 2) 
E3( 3) 


C 

C- 

c 

c- 

c 


c 

c- 

c 

c 


c 

c- 

c 

c 

c 

n 

c 

c 

c 

c 

C 

c 


GO TO 1000 

MATED DUCT SURFACE 

600 I HOLE - 1 

EPS1 » ETA(l) 

EPS3 - ETA( 3 ) 

— SUB-SURFACE CORNER POINTS 

DO 610 J-1,6 

R1(J) « PL(J,4,2) 

R2(J) - PL( J ,5,2) 

R3( J) - PL(J#5,3) 

610 R4(J) - PL(J,4,3) 

— SUB-SURFACE EDGES 

B ■ BETA1 + EPS1* ( BETA2 

CALL HOLE (B, El, VECTOR) 

B - BETA2 + EPS 3* ( BETA3 

CALL HOLE (B,E 2, VECTOR) 

B » BETA4 + EPS1* (BETA3 

CALL HOLE (B,E 3, VECTOR) 

B « BETA1 + EPS3MBETA4 

CALL HOLE (B,E 4, VECTOR) 

INTERPOLATION 


BETA1 ) 

BETA2 - 2.*PI) 
BETA4 ) 

BETA1 ) 


MAP ■ 6 
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— POSITION 

DO 620 J»l,3 

620 RAM(J) ■ (1. - EPS1)*(E4(J) - (1. - EPS3)*R1(J) - EPS3*R4(J)) 

1 + EPS1* (E2( J) - (1. - EPS3)*R2(J) - EPS3*R3(J)) + E1(J) 

2 + EPS3*(E3(J) - El ( J ) ) 


RADIUS AND TANGENT 
— AXIAL DISTANCE 

RAX - RAM(l) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


-RADIUS 

IF( EPS1.EQ.0. .OR. EPS 1 . EQ. 1 . .OR.EPS3.EQ.O. .OR. EPS3. EQ. 1 . ) THEN 

RAD » DFNR + BE*SQRT(1.0 - ((RAX - DFND)/AE)**2) 

RAD » SORT ( RAM ( 2 ) **2 + RAM(3)**2) 

ELSE 

I - 0 

THETA ■ AT AN ( RAM( 2) / RAM( 3) ) 

DO 640 J»1 , 36 

CALL HOLE ( BETA (J), El, VECTOR) 

THETA1 » ATAN( El(2) / El(3) ) 

CALL HOLE(BETA( J+l ) ,E2, VECTOR) 

THETA2 - AT AN ( E2(2) / E2(3) ) 

MATCH ■ 0 

I F ( THETA . GT . THETA1 .AND. THETA. LT.THETA2) MATCH » 1 
IF( THETA. GT.THETA2 .AND. THETA. LT.THETA1 ) MATCH - 2 

IF ( MATCH. GT.O) THEN 

I ■ I + 1 

I F ( MATCH . EQ • 1 ) THEN 

RATIO • (THETA - THETA1) / (THETA2 - THETA1) 

B » RATIO* (BETA( J+l) - BETA(J)) + BETA(J) 

ELSE 

RATIO - (THETA - THETA2) / (THETA1 - THETA2 ) 

B * RATIO* (BETA(J) - BETA (J+l ) ) + BETA(J+1) 
END IF 

IF(I.EQ.l) THEN 

CALL HOLE (B,E3, VECTOR) 
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C 

c 

c 

c 

c 

c 

c 

c 

640 CONTINUE 
650 CONTINUE 
C 

IF(E3( 1) 
C 

C 

C 

C 

c 

c 


RADI - S0RT( E3( 2)**2 + E3(3)**2) 

GO TO 640 
ELSE 

CALL HOLE ( B r E 4 r VECTOR ) 

RAD 2 ■ SQRT( E4(2)**2 ♦ E4(3)**2) 

GO TO 650 
END IF 

END IF 

LT.E4 ( I ) ) THEN 

RATIO ■ (RAM(l) - E3( 1 ) ) / (E4(l) - E3(l)) 

RAD ■ RATIO* (RAD2 - RADI) + RADI 
ELSE 

RATIO - (RAM(l) - E4( 1) ) / (E3{1) - E4(l)) 

RAD - RATIO* (RADI - RAD2) + RAD2 
END IF 

END IF 


— TANGENT 


DRDX - -( ( BE/AE ) **2 ) * ( RAX - DFND)/(RAD - DFNR) 
TH ■ AT AN ( DRDX ) 

DSDN (1,1) « COS(TH) 

DSDN (2,1) - SIN(TH) 

DSDN (3,1) - 0.0 


OUTPUT POSITION AND NORMAL 


-AXIAL COMPONENT OF TANGENT 
DA ■ DSDN (1,1) 

-NORMAL COMPONENT OF TANGENT 

DN - DSDN( 2,1) 

-ANGULAR LOCATION 

RADX « SORT(RAM(2)**2 + RAM(3)**2) 
ANG ■ ASIN(RAM(3)/RADX) 

CANG - COS (ANG) 

SANG * SIN (ANG) 
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C 

C TANGENT 1 (AXIAL DIRECTION) 

C 

DSDN (1,1) * DA 
DSDN (2*1) - DN*CANG 
DSDN (3,1) « DN*SANG 
C 

C TANGENT 2 (CIRCUMFERENTIAL DIRECTION) 

C 

DSDN (1,2) - 0.0 

DSDN( 2,2) » -SANG 
DSDN( 3,2) - CANG 
C 

C NORMAL 


C 

c 

c- 

c 


c 

c 


c 


CALL CROSS ( DSDN (1,1), DSDN (1,2), SNORMAL ( 1 , ISIDE ) , 46 ) 
-POSITION 


SIDE(1, ISIDE) 
SIDE(2, ISIDE) 
SIDE(3, ISIDE) 


RAX 

RAD* CANG 
RAD* SANG 


IF(EPS3.EQ.1.0) THEN 

SIDE(1, ISIDE) - E3d ) 
SIDE(2, ISIDE) » E3( 2) 
SIDE(3, ISIDE) » E3( 3) 


END IF 


C 

c- 

c 

c- 

c 

c- 

c 


c 

c 

c 

c 


c 

c 

c 

c 


GO TO 1000 

DUCT SURFACE NEAR INLET 


700 IF( ISIDE.LT. 4) THEN 
—SURFACE 1 AND 3 

IEDG2 - LINE4 

EPS1 ■ ETAd) 
EPS2 » ETAd) 

ELSE 

— SURFACE 5 AND 6 

IEDG2 - LINE3 

EPS1 » ETAd) 
EPS 2 ■ ETAd) 

END IF 

WELD REGION 

PI2 - 1.57079633 


MAP - 7 


C 

B « COEFS(l, ISIDE) ♦ EPSl*(COEFS(2,ISIDE) - COEFSd, ISIDE) ) 


EPS 3 - (COEFSd, ISIDE) - 1.0)/(NMBRNDS( 2) - 1.0) 
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C 


RATIO « EPS2/EPS3 
IF ( RATIO. GT. 1.0) RATIO « 1.0 


—WELD RADIUS (RDOC) i DIFFERENCE BETWEEN HOLE AND DUCT RADIUS ( DRHO ) 


CALL DELRAD ( B , RHO , DRHO , RHOD , RDOC ) 

— CALCULATE AXIAL DISTANCE (CP) 

CALL CAXIS(B,RHO,CP) 

Y « RDOC - DRHO 

AAD » 0.5*( (ABOT + ATOP) + (ATOP - ABOT)*COS(B - PI2)) 
X « SQRT(DRHO*( 2.0* RDOC - DRHO)) 

— WELD ANGLE 

TH2 » ATAN(X/Y) 

TH1 ■ PI2 - TH2 


C- 

c- 

c~ 

c~ 


XX » AAD - X 
TH3 ■ ATAN ( XX/ RDOC ) 

ANGL * RATIO* (TH2 + TH3) + TH1 
— POSITION AND TANGENT 1 (AXIAL DIRECTION) 
IF( ANGL.GE.PI2) THEN 

— DUCT 

-AXIAL DISTANCE 
-RADIUS 


-COORDINATES 

-TANGENT 


PSI » (ANGL - PI2J/TH3 
CS ■ CP ♦ X ♦ PSI*XX 
RS * RHOD 

CALL DUCT(B,RS,CS,SIDE(l r ISIDE)) 


DSDN(1,1) * DUl(l) 
DSDN (2,1) - DU1(2) 
DSDN (3,1) ■ DU1(3) 


C WELD- 


ELSE 


C 

C AXIAL DISTANCE 

C RADIUS 

C COORDINATES 

C TANGENT 

C 


XS - RDOC*COS ( ANGL ) 

YS « RDOC*SIN( ANGL) 

CS » CP + X - XS 

RS • RHOD + RDOC - YS 

CALL DUCT(B,RS,CS,SIDE(1,ISIDE)) 

DXDN • TAN (ANGL) 

CB • COS(B) 

SB • SIN(B) 
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C 

CALL VADD(CB, DU2,SB,DU3, VECTOR, DSDN(1, 1) ) 

C 

ONE - 1.0 

C 

CALL VADD ( DXDN , DU 1 , -ONE , VECTOR r C , DSDN (1,1)) 

C 

END IF 

C 

C DUCT REGION 

C 

IF( RATIO. EQ* 1.0) THEN 

C 

C POSITION 

PSI « (EPS 2 - EPS3 )/( 1 • - EPS 3 ) 

PSI1 « 1.0 - PSI 

C 

CALL VADD ( PS II , SIDE( 1 , IS IDE ) , 

& PSI, EDGE ( 1 , 1 EDG 2 ) , S IDE ( 1 , IS IDE ) , VECTOR ) 

C 

c TANGENT 1 (AXIAL DIRECTION) ; INPUT 

C 

DSDN (1,1) • DU 1 ( 1 ) 

DSDN (2,1) « DUK 2) 

DSDN (3,1) - DUK 3) 

C 

END IF 

C 

C TANGENT 2 (CIRCUMFERENTIAL DIRECTION) 

C 

CB « COS(B) 

SB * SIN ( B ) 

C 

CALL VADD ( CB , DU 2 , SB , DU3 , VECTOR , R1 ) 

C 

CALL CROSS ( DSDN (1,1) , VECTOR, DSDN ( 1 , 2 ) , 47 ) 

C 

C NORMAL 

C 

CALL CROSS ( DSDN (1,1), DSDN (1,2) , SNORMAL( 1 , IS IDE) , 48 ) 

C 



C DIRECT SURFACE NORMAL INTO FLOW DOMAIN 



1000 DO 1010 1*1 , 3 

1010 SNORMAL( I , IS IDE) * SNORMAL( I , ISIDE) # SIGNS( ISIDE ) 

C 

RETURN 

C 

C FORMAT STATEMENTS 

C 

1100 FORMAT ( 1H1 , 10X, 20A4 , 13X, 8H SECTION, 12 , 3H OF,I3,9H FOR ZONE, 13) 

C 

END 

C 

c * ****************************************************************** 
C ***** 0UTPUT * ********** ********************************* ************ 
c ******************************************************************* 

C 

SUBROUTINE B LKOUT ( NUN IT, NODSTOR ) 
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C 

C WRITES THE FORMATTED BLOCKED GEOMETRY FILE (NUNIT) 

‘ 

COMMON /INITA/ I ZINDEX,MAPTEN, INCHES 

COMMON /IOCOUNT/ IREWINDt 40 ) ,NREAD( 40 ) ,NWRITE( 40 ) 

COMMON /MARCHS/ MARCH , INDEX( 3 ) 

COMMON /OUT/ NODE (3,4000) 

COMMON /ZONING/ I ZONE , ISECT , NZ INDEX , NMBRNDS ( 3 ) 

C 

IPLANE « NWRITE( NUNIT) + 1 
C 

WRITE(NUNIT, 1000) NODSTOR, IPLANE, 

& NMBRNDS ( 1 ) , NMBRNDS ( 2 ) , NMBRNDS ( 3 ) , MARCH 

WRITE( NUNIT, 1010) ( NODE ( 1 , I ) , I s 1 , NODSTOR) 
WRITE( NUNIT, 1010 ) (NODE( 2 , I ) , 1-1 , NODSTOR) 
WRITE ( NUNIT , 1010 ) (NODE(3, I) , I»l, NODSTOR) 
C 

10 CONTINUE 
C 

40 NWRITE( NUNIT) ■ NWRITE ( NUNIT) + 1 
C 

RETURN 

C 

C FORMAT STATEMENTS 

C 

1000 FORMAT( 2415) 

1010 FORMAT( 6E22. 14 ) 

1020 FORMAT( 2216 ) 

C 

END 


C 

c ****************************************************< 

**************************** 

c ****************************************************1 

c 

SUBROUTINE CROSS(A,B,C,N) 

C 

C C * CROSS PRODUCT OF A AND B (UNIT VECTOR) 

C 

COMMON /COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
COMMON /HEADER/ ITITLE(20) , LINE 
COMMON /INITA/ IZINDEX,MAPTEN, INCHES 
COMMON /UNITS/ NU5,NU6,NU20 

COMMON /ZONING/ I ZONE, ISECT, NZ INDEX, NMBRNDS ( 3 ) 

C 

DIMENSION A( 3) ,B( 3) ,C( 3 ) 

C 

C CROSS PRODUCT 

C 

C(l) - A(2)*B<3) - A(3)*B(2) 

C ( 2 ) - A(3)*B(1) - A(1) # B(3) 

C( 3 ) - A( 1 )*B( 2) - A(2)*B(1) 

C 

C MAGNITUDE 

C 

CALL VMAG(C,CMAG) 

C 

IF(CMAG.GT.O.O) THEN 

C 

C NORMALIZE 
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C 


C 


C 


C(l) - C ( 1 ) /CMAG 
C ( 2 ) - C ( 2 ) / CMAG 
C( 3) - C { 3 ) /CMAG 

ELSE 

C(l) * 0.0 
C(2) * 0.0 
C( 3) » 0.0 

END IF 

RETURN 


—FORMAT STATEMENTS 

000 FORMAT ( 1H1 , 10X, 20A4 , 13X, 8H SECTION, 12 , 3H OF,I3,9H FOR ZONE, 13) 

010 FORMAT ( 9H LOCATION, 13, 36H: CROSS PRODUCT EQUALS ZERO FOR NODE, 16) 

END 


GRID SPACING 


SUBROUT I N E ETABC ( MARCH , INDEX , NODE ) 


THIS ROUTINE SEPERATES THE BOUNDARY CONDITIONS AND 
CALCULATES THE VALUE OF ETA. 


COMMON /INITA/ I Z INDEX, MAPTEN, INCHES 
COMMON /INITC/ PI ,RADDEG 

COMMON /INPUTA/ EDGE( 3,12) , POINT( 3, 8 ) ,SIDE( 3,6) 

COMMON /INPUTBC/ INODEBC( 3 ) , ISIDE( 3 ) 

COMMON /SPACING/ ISTRTCH(3) , STRETCH ( 3 ) , ETAS (3,200) , ETA( 3 ) , DETA( 3 ) 
COMMON /ZONING/ IZONE, ISECT,NZ INDEX, NMBRNDS( 3 ) 


FIRST NODE 


IF(NODE.EQ.l) THEN 

ETA( INDEX) » 0.0 


— DETERMINE SIDE 


ISIDE( INDEX) « 10 + (INDEX - 6)*INDEX 

STORE SPACING 


ETAS ( INDEX, 1) - 0.0 

RETURN 

END IF 


C 

C 


LAST NODE 

IF ( NODE. EQ.NMBRNDS( INDEX)) THEN 

ETA( INDEX) >1.0 
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— DETERMINE SIDE 

ISIDE( INDEX) - 9 + (INDEX - 7)*INDEX/2 

— STORE SPACING 

ETAS ( INDEX , NODE ) « 1.0 

RETURN 

END IF 

INTERIOR NODES 

IN0DEBC( INDEX) » 9 

IF( ISECT.GT. 1 .AND. I NDEX . NE . MARCH ) GO TO 310 
— CALCULATE ETA 

ISTR - ISTRTCHt INDEX) + 1 

GO TO (100,110,120,130,140,150,160,170,180,190,200) ISTR 


— EQUAL SPACING (0) 

100 ETA( INDEX) « ETA( INDEX) + DETA( INDEX) 

GO TO 300 

— INPUT ETA SPACING (1) 

110 ETA (INDEX) - ETAS ( INDEX, NODE) 

GO TO 300 

— DECREASING SPACING; INPUT STRETCHING FACTOR (2) 


120 RATIO - REAL ( NODE - 1 ) /REAL (NMBRNDS( INDEX) - 1) 

XI » RATIO*STRETCH( INDEX )/2.0 
ETA1 * (EXP(Xl) - EXP ( -XI ) ) / ( EXP ( XI ) + EXP(-Xl)) 

X2 • STRETCH (I NDEX)/ 2.0 

ETA 2 - ( EXP ( X2 ) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 

ETA( INDEX) » ETA1/ETA2 
GO TO 300 

INCREASING SPACING; INPUT STRETCHING FACTOR (3) 

130 RATIO « REAL(NMBRNDS( INDEX) - NODE ) /REAL (NMBRNDS( INDEX) - 1) 

XI ■ RATIO*STRETCH( INDEX) /2.0 
ETA1 - (EXP(Xl) - EXP(-Xl) )/(BXP(Xl) + EXP(-Xl)) 

X2 « STRETCH ( INDEX )/2.0 

ETA 2 • ( EXP ( X2 ) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 
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ETA (INDEX) ■ 1.0 - ETA1/ETA2 

C 

GO TO 300 

— DOUBLE STRETCHING; INPUT STRETCHING FACTOR (4) 

140 XI - ( STRETCH ( INDEX)/2.0)/2.0 

ETAMID « (EXP (XI) - EXP ( -XI ) )/( EXP( XI ) + EXP(-Xl)) 

RATIO - REAL (NODE - 1 ) /REAL (NMBRNDS( INDEX) - 1) 

X2 « (RATIO - 0 . 5 ) ‘STRETCH ( INDEX)/2.0 
ETAMAX ■ (EXP(X2) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 

ETA( INDEX) - (ETAMAX + ETAMID) /( 2. 0*ETAMID) 

GO TO 300 

— DECREASING SPACING; INPUT MINIMUM SPACING (5) 

150 ARGI - STRETCH (INDEX) ‘REAL (NODE - 1) 

EXPI - EXP (ARGI ) 

EXP I I - 1 . 0/EXP I 

TAN HI « (EXPI - EXPII )/( EXPI + EXPII) 

ARGN - STRETCH ( INDEX) *REAL(NMBRNDS( INDEX) - 1) 

EXPN ■ EXP (ARGN) 

EXPNI » 1.0/EXPN 

TANHN - (EXPN - EXPNI) /(EXPN + EXPNI) 

ETA( INDEX) « TANHI/TANHN 
GO TO 300 

— INCREASING SPACING; INPUT MINIMUM SPACING (6) 

160 ARGI » STRETCH ( INDEX) ‘REAL ( NMBRNDS ( INDEX) - NODE) 

EXPI * EXP (ARGI ) 

EXPII « 1 . 0/EXP I 

TANHI » (EXPI * EXPI I )/( EXPI + EXPII) 

ARGN - STRETCH( INDEX) ‘REAL (NMBRNDS (INDEX) - 1) 

EXPN - EXP (ARGN) 

EXPNI - 1.0/EXPN 

TANHN « (EXPN - EXPNI)/ (EXPN ♦ EXPNI) 

ETA( INDEX) - 1.0 - TANHI/TANHN 
GO TO 300 

— DOUBLE STRETCHING; INPUT MINIMUM SPACING (7) 

170 ARGI • STRETCH ( INDEX) *REAL(2‘NODE - NMBRNDS ( INDEX) - 1) 

EXPI * EXP (ARGI) 

EXPII » 1.0/EXPI 

TANHI ■ (EXPI - EXPII )/ ( EXPI + EXPII) 

ARGN * STRETCH( INDEX)*REAL(NMBRNDS( INDEX) - 1) 

EXPN - EXP (ARGN) 

EXPNI « 1.0/EXPN 
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TANHN • (EXPN - EXPNI)/(EXPN + EXPNI) 

ETA( INDEX) - 0.5* (1.0 + TANHI/TANHN) 

GO TO 300 

—DECREASING SPACING; INPUT STRETCHING FACTOR 

180 PIDN - PI/ ( STRETCH ( INDEX) *NMBRNDS ( INDEX) ) 

ETA( INDEX) - 1.0 - TAN ( PIDN* ( NMBRNDS ( INDEX) - NODE)) 
6 /TAN (PIDN* (NMBRNDS (INDEX) -1)) 

GO TO 300 

— INCREASING SPACING; INPUT STRETCHING FACTOR 

190 PIDN ■ PI/ ( STRETCH ( INDEX) *NMBRNDS( INDEX) ) 

ETA( INDEX) - TAN ( PIDN* ( NODE - 1)) 

6 /TAN ( PIDN* ( NMBRNDS (INDEX) - 1)) 

GO TO 300 

USER INPUT STRETCHING FUNCTION 

200 CONTINUE 

STORE SPACING 

300 ETAS (INDEX, NODE) » ETA (INDEX) 

310 ETA (INDEX) * ETAS ( INDEX, NODE) 

RETURN 

END 

r **********************•********************************< 
‘OUTPUT****** **************************************** 


C 

C 

C 

C- 


SUBROUTINE IOPACK(NUNIT) 

GENERAL PURPOSE FORTRAN I/O PACKAGE FOR UNITS 1 -> 40 

COMMON /IOCOUNT/ IREWIND( 40 ) ,NREAD( 40 ) ,NWRITE( 40 ) 
COMMON /UNITS/ NU5,NU6,NU20 

DATA I REWIND , NREAD , NWRITE /120*0/ 

ENTRY RWIND: REWIND FILE ON NUNIT 

ENTRY RWIND( NUNIT) 

IREWIND( NUNIT) ■ 1 

REWIND NUNIT 

RETURN 


(8) 


(9) 


(10) 
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ENTRY STATUS: PRINT STATUS OP I/O OPERATIONS ON ALL UNITS 


ENTRY STATUS (NUNIT) 

WRITE (NU6, 1000) 

DO 10 NU-1,40 

10 IF( IREWIND(NU) . EQ.l) WRITE( NU6 , 1010 ) NU , NREAD ( NU ) , NWRITE ( NU ) 
RETURN 

-FORMAT STATEMENTS 

00 FORMAT (//42X,38H STATUS OF I/O OPERATIONS ON ALL UNITS 

1 //42X, 32H UNIT NO. OF NO. OF • 

2 /42X, 32H NUMBER READS WRITES / ) 

1010 FORMAT ( 4 2X , 3 1 1 0 ) 

END 


****OUTPUT 


SUBROUTINE OUTPUT ( NUN IT, NODSTOR) 
PRINTOUT AND STORE DATA 


COMMON /COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
COMMON /HEADER/ ITITLE( 20 ) , LINE 
COMMON /INITA/ 12 INDEX, MAPTEN , INCHES 
COMMON /OUT/ NODE (3,4000) 

COMMON /UNITS/ NU5,NU6,NU20 

COMMON /ZONING/ IZONE,ISECT,NZINDEX,NMBRNDS(3) 


-TOTAL NUMBER OF PLANES 

NPLANE - NPLANE + 1 


C 

c 

c 

c 

c 


PRINT NODAL INFORMATION 

IPRINT - 1 

NODETOT - NODSTOR + NODETOT 
NODESAV » NODESAV - NODSTOR 

PRINT TOTAL NUMBER OF POINTS STORED 

LINE - LINE + 3 

IF( LINE.GE. 60) THEN 

WRITE (NU6, 1000) ITITLE,IZONE 
WRITE (NU6, 1030) 

LINE - 7 

END IF 
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WRITE (NU6,1060) NODSTOR , NPLANE , NUN IT , NODETOT 


STORE OUTPUT 


200 CALL BLKOUT(NUNIT, NODSTOR) 
RETURN 


FORMAT STATEMENTS 


L000 

L020 

L030 

L040 

1050 


1060 


FORMAT ( 1H1 , 10X, 20A4 , 13X, 8H SECTION , 12 , 3H OF,I3,9H FOR ZONE, 13) 
FORMAT (1X,I6,2(1X,F13.7),5X,F7.2»6X»F7.2,5X,I2) 

FORMAT ( / 44H NODE X Y Z ) 

FORMAT ( 1X,I6,3(1X,F13.7),2X,2(2X,F7.2),F11.3,3X,2(2X,F7.2),5X,I2) 
FORMAT ( / 5X, 15 , 18H POINTS FROM PLANE, 14, 15H STORED ON UNIT, 13, 

23H: TOTAL POINTS STORED «,I6 / ) 

FORMAT (/ 10X, 15, 18H POINTS FROM PLANE, 14, 15H STORED ON UNIT, 13, 
23H: TOTAL POINTS STORED «,I6 / ) 


END 


*****OUTPUT 


SUBROUTINE PICTURE( IDRAW) 


THIS ROUTINE DESCRIBES THE NOMENCLATURE 


COMMON /UNITS/ NU5,NU6,NU20 

WRITE ( NU6 , 300 ) 

WRITE(NU6,310) 

WRITE ( NU6 , 320 ) 

WRITE(NU6,330) 

RETURN 

FORMAT STATEMENTS 


300 FORMAT! 40X,37H 

1 /40X, 19H 

2 /40X, 19H 

3 /40X, 34H 

4 /40X, 32H 


5 / 40X, 50H 

6 /40X, 42H 

7 /40X, 42H 

8 /40X, 50H 

310 FORMAT ( 40X,50H 


1 

2 

3 

4 

5 

6 

7 

8 
9 


/40X, 50H 
/40X, 50H 
/40X, 50H 
/40X, 50H 
/40X, 50H 
/40X, 46H 
/ 40X, 50H 
/40X,42H 
/40X, 50H 


NOMENCLATURE 

IE 

IT 

IA SURFACE 4 

12 (TOP) 



POINT 4 

I 



POINT 
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0— 
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320 FORMAT ( 40X,51H 

1 /40X, 42H 

2 / 40X, 41H 

3 /40X, 40H 

4 /40X, 39H 

5 /40X, 38H 

6 /40X, 44H 

7 /40X, 10H 

8 /40X.26H 

9 /40X, 25H 
330 FORMAT ( 40X, 8H 

6 /40X, 7H 


C POINT 1/5 EDGE 1 

11 / POINT 2 

Ell /E 

10 /6 

2 I /G 

I /E 

5 I /D 

I /G 

I /E 

I /D 

0/ 

O/E 

POINT 5 EDGE 9 

POINT 6 


/ 

/3 SURFACE 2 

/A (BOTTOM) ) 

/T 

/E ) 


123456789012345678901234567890123456789012345678901 


END 


•‘•‘UTILITY 


SUBROUTINE VADD(CA,A,CB,B,C,UC) 


VADD COMPUTES C,THE SUM OF VECTORS CA*A AND CB*B, WHERE CA AND 
CB ARE SCALARS. UC IS A UNIT VECTOR DIRECTED ALONG C. 


DIMENSION A(3) ,B(3) ,C(3) ,UC(3) 
SUM * 0.0 
DO 10 1*1,3 

C(I) = CA*A( I ) + CB*B( I ) 

10 SUM * SUM + C(I)*C(I) 

CMAG * SORT (SUM) 


RMAG * 0.0 

I F ( CMAG . GT .0.0) RMAG - l.O/CMAG 

UC(1) * C ( 1 ) ‘ RMAG 
UC(2) - C ( 2 ) * RMAG 
UC( 3) - C ( 3 ) *RMAG 

RETURN 

END 

♦a********************************************************************* 
**** ‘UTILITY** ********************************************************* 


SUBROUTINE VDOT(A,B,C) 

VDOT COMPUTES C,THE DOT PRODUCT OF VECTORS A AND B. 


DIMENSION A(3) ,8(3) 
C 

C « 0.0 

c 
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DO 10 1-1,3 
10 C * C + A( I ) *B{ I ) 

RETURN 

END 


**** UTILIT y*»********************»******»******************** 
ft************************************************************ 1 

SUBROUTINE VMAG( VECTOR, VECMAG) 

VMAG DETERMINES THE MAGNITUDE OF A VECTOR 

DIMENSION VECTOR ( 3 ) 

VECMAG « SORT ( VECTOR ( 1 ) * * 2 + VECTOR(2)**2 + VECT0R( 3 ) **2 ) 

IF(VECMAG.LT. 0.0000001) VECMAG ■ 0.0 

RETURN 

END 


****HGM 2******************** 
****** *T* *************** ***** 

SUBROUTINE CAXIS(B,RHO,C) 


CALCULATE AXIAL DISTANCE CORRESPONDING TO A POINT ON THE HOLE. 

B ANGULAR LOCATION OF POINT ON THE HOLE 

C AXIAL DISTANCE 

RHO HOLE RADIUS 


COMMON /DFN1/ DU1 ( 3 ) , DU2 ( 3 ) , DU3 ( 3 ) ,DFNB,DFND, DFNF 
COMMON /DFN2/ AE , BE , DFNR, ZSTAR, AGL,BETA1 , BETA2 , BETA3 ,BETA4 
COMMON /INITC/ PI ,RADDEG 

—TOLERANCE FOR NEWTON-RAPHSON ITERATION 

FEPS « 1.0E—07 

— ANGLE OF DUCT AXIS 

SA - SIN(AGL) 

CA - COS (AGL) 

ANGLULAR LOCATION OF A POINT ON THE HOLE 

SB « SIN(B) 

CB - COS(B) 

AXIAL DISTANCE ALONG LOWER EDGE OF DUCT 

DELTA - 0 . 0/RADDEG 

C 

IF( (B.GE. (BETA4 - DELTA)) .AND. (B.LE.BETA3) ) THEN 
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C 

C 

C 


C ■ ( ZSTAR + RHO*SA*CB)/CA 
RETURN 
END IF 


— DISTANCE ALONG MAJOR AXIS OF BOWL ELLIPSE 

XT - DFNF - DFND + RHO*CB 
I F ( ABS ( XT ) . GT . AE ) XT - AE 

— VERTICAL DISTANCE TO HOLE 

YT • DFNB + RHO*SB 

— DISTANCE ALONG MINOR AXIS OF BOWL ELLIPSE 

DR » SQRT(AE**2 - XT**2)*BE/AE 

— INITIAL DUCT AXIAL DISTANCE 

C » SORT ( ( DFNR + DR)**2 - YT**2) 

— NEWTON-RAPHSON ITERATION 

10 XTT - DFNF - DFND + C*SA + RHO*CA*CB 
YTT - YT 

ZTT - C*CA - RHO*SA*CB 
— FUNCTION 

DR - SQRT(AE**2 - XTT**2 ) *BE/AE 

FC - (DFNR + DR) **2 - YTT**2 - ZTT**2 

— DERIVATIVE 

DFDC1 - XTT*SA* (DFNR + DR) *BE/AE 
DFDC2 - DR*AE/BE 
DFDC3 « ZTT*CA 

DFDC ■ -2 . * { DFDC1/DFDC2 + DFDC3 ) 

— DUCT AXIAL DISTANCE 

C - C - FC/DFDC 

IF(ABS(FC) .GT.FEPS) GO TO 10 

RETURN 

END 

►****HGM 2************************************** 
>**•*•*•▼**«•****** ***************************** 

SUBROUTINE DCTDAT 

DATA DESCRIBING HOLE AND WELD RADIUS 
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COMMON /DFN3/ AD, BD,BETA( 37 ) , DELRHO( 37 ) ,RADOC ( 37 ) 

C 

DIMENSION DCT1 ( 37 ) , DCT2( 37 ) ,DCT3 ( 37 ) 

C 

DATA RADDEG/57. 29577951/ 

—ANGULAR INCREMENT AROUND HOLE 


—DIFFERENCE IN RADIUS BETWEEN DUCT AND HOLE 


— RADIUS OF CURVATURE OF FAIRING 


— READ DUCT DATA 

READ( 5,1000) (DCTK I), 1-1,37) 

READ( 5 ,1000 ) (DCT2(I) ,1-1,37) 

READ( 5 , 1000 ) (DCT3( I) ,1-1,37) 

— CONVERT DEGREES TO RADIANS AND INCHES TO FEET 

DO 10 I - 1,37 

BETA(I) » DCTK I) / RADDEG 
DELRHO(I) - DCT2 ( I ) 

10 RADOC(I) - DCT3( I) 

— FORMAT STATEMENT 

.000 FORMAT ( (8F10.3) ) 

RETURN 

END 


****HGM 2 
****** *T * 


SUBROUTINE DELRAD ( B , RHO , DRHO , RHOD , RDOC ) 


INTERPOLATES HOLE RADIUS, WELD RADIUS, DUCT RADIUS, AND DIFFERENCE 
BETWEEN HOLE AND DUCT RADIUS FOR POINTS BETWEEN ANGULAR INCREMENTS. 

B ANGULAR LOCATION OF POINT ON HOLE 

RHO HOLE RADIUS 

RHOD DUCT RADIUS 

RDOC WELD RADIUS 

DRHO DIFFERENCE BETWEEN HOLE AND DUCT RADIUS 


COMMON /DFN3/ AD,BD,BETA( 37) ,DELRHO( 37 ) ,RADOC( 37 ) 
— ANGULAR LOCATION OF A POINT ON THE DUCT 
BB — B 

IF(B.LT.O.O) BB - 6.283185308 + B 


C-60 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 


non noon oonnn o n n non n non non non noon n n non 


LMSC-HEC TR F225968 


BRACKET ANGULAR INCREMENT 

DO 10 I - 1,36 

DBK * B£TA( 1+1 ) - BETA(I) 

DB - BB - BETA(I) 

TF(DB.LE.DBK) THEN 

INTERPOLATE 

SEGMENT RATIO 

PS I - DB/DBK 

DIFFERENCE IN RADIUS BETWEEN DUCT AND HOLE 

DRHO - DELRHO(I)*(l. - PSI) + DELRHO( 1+1 ) *PSI 

RADIUS OF WELD 

RDOC - RADOC(I)*(l. - PSI) + RADOC(I+l)*PSI 

RADIUS OF DUCT ELLIPSE 

CBSQ - COS(B)**2 
SBSQ - SIN(B)**2 

RADX p SORT ( SBSQ*AD*AD ♦ CBS0*BD*BD) 

RHOD » AD*BD/RADX 

HOLE RADIUS 

RHO - RHOD + DRHO 
RETURN 
END IF 

10 CONTINUE 

RETURN 

END 

*****HGM 2 ****************************************************** 

»****««*T** ******* *•**#**•**• ********** •****••** **************** 

SUBROUTINE D£RIV( B , RHO , C , DRDB , DCDB ) 

COMPUTE DRDB, DCDB 

COMMON /DFN1/ DU1 (3) ,DU2( 3) ,DU3 (3) ,DFNB,DFND,DFNF 
COMMON /DFN2/ AE , BE , DFNR , ZSTAR , AGL , BETA1 , BETA2 , BETA3 , BETA4 
COMMON /DFN3/ AD, BD, BETA( 37 ) ,DELRHO( 37 ) ,RADOC( 37 ) 

POINT ANGLE 

CB « COS(B) 

SB * SIN(B) 
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EX1 « BD**2 - AD* *2 

EX2 - ( (AD*SB)**2 + (BD*CB)**2)**1.5 
DRDDB - EX1*AD*BD*SB*CB/EX2 

BRACKET BETA 

DO 10 I « 1,36 

DBI « BETA( 1+1 ) - BETA(I) 

DB « B - BETA(I) 

IB » I 

IF(DB.LE.DBI) GO TO 20 

10 CONTINUE 
20 CONTINUE 

DDRO » (DELRHO( IB+1) - DELRHOC IB) )/(BETA( IB+1 ) - BETA(IB)) 

DRDB » DRDDB + DDRO 

FOR BETA BETWEEN BETA4 AND BETA3 

I F { (B.GE.BETA4) .AND. (B.LE.BETA3) ) THEN 

DCDB ■ ( DRDB*CB-RHO*SB) *TAN ( AGL) 

RETURN 

END IF 


AXIS ANGLE 

CA - COS (AGL) 

SA - SIN(AGL) 

XMD ■ DFNF - DFND ♦ C*SA + RHO*CA*CB 
EXR1 « BE/AE 

EXR2 ■ SQRT( AE**2 - XMD**2) 

EXR - (DFNR + EXR1*EXR2 ) *XMD*EXR1/EXR2 

EXC1 - C*CA - RHO*SA*CB 
EXC2 ■ RHO*CB + DRDB*SB 
EXC3 - DRDB*CB - RHO*SB 
EXC4 - DFNB + RHO*SB 

EXCN - ( EXC1*SA - EXR*CA) *EXC3 - EXC2*EXC4 
EXCD » EXR*SA ♦ EXC1*CA 
DCDB » EXCN/EXCD 

RETURN 

END 


*****HGM 2*************************************** 
********T ************* ****** ***** **************** 

SUBROUTINE DEXIT(B, POINT, DE) 

COMPUTE POINT COORDINATES ON DUCT EXIT PLANE 
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AGL ANGLE OF DUCT AXIS 

B ANGLE OF POINT 

AD ELLIPSE MINOR AXIS (HORIZONTAL) 

BD ELLIPSE MAJOR AXIS (VERTICAL) 

CO AXIAL DISTANCE TO POINT ON DUCT EXIT NORMAL IN X Z PLANE 

CE AXIAL DISTANCE TO POINT ON AXIS OF DUCT EXIT 

DFNB VERTICAL DISTANCE TO AXIS ORIGIN 
DFNF AXIAL DISTANCE TO AXIS ORIGIN 


COMMON /DFN1/ DU1 ( 3 ) , DU2 ( 3 ) ,DU3 ( 3 ) , DFNB ,DFND, DFNF 
COMMON /DFN2/ AE , BE , DFNR , ZSTAR , AGL , BETA1 , BETA2 , BETA3 , BETA4 
COMMON /DFN3/ AD, BD, BETA< 37 ) ,DELRHO( 37 ) ,RADOC( 37 ) 

COMMON /DFN5/ CE,CQ,ABOT,ATOP 

DIMENSION POINT ( 6 ) 

— AXIS ANGLE 

COA - COS (AGL) 

SIA « SIN(AGL) 

POINT ANGLE 

COB - COS(B) 

SIB ■ SIN(B) 

— RADIUS OF ELLIPTIC DUCT 

RADX - SQRT( (AD*SIB)**2 + (BD*COB)**2) 

RHOD ■ AD*BD/RADX 

—AXIAL DISTANCE 

C ■ CE + RHODES IB*DFNB/( CO - CE) 

— COORDINATES 

POINT(l) « DFNF + C*SIA + RHOD*COA*COB 
POINT ( 2 ) * DFNB + RHOD*SIB 

POINT( 3 ) « C*COA - RHOD*SIA*COB 

RETURN 

END 


***** HGH 2 
******* 


SUBROUTINE HGMIN 


INPUT FOR TWO DUCT HGM WITH TURN AROUND DUCT 


COMMON /COEFF/ 
COMMON /INITA/ 
COMMON /INITC/ 

COMMON /INPUTA/ 
COMMON /MAPING/ 


COEFE( 8,10 ,12) ,COEPS( 8 ,6) ,NMBRSEG( 12) 

I Z INDEX , MAPTEN , INCHES 

PI,RADDEG 

EDGE( 3,12) , POINT ( 3,8 ) , SIDE ( 3,6) 
MAPSIDE( 6 ) , MAPS EG ( 10,12) 
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COMMON /MAXIMUM/ ETAMAX( 10 , 12 ) , SEGMAX( 3 , 10 , 12 ) 

COMMON /ZONING/ IZONE, ISECT,NZINDEX,NMBRNDS( 3) 

C 

COMMON /DFN1/ DUI ( 3 ) , DU2 ( 3 ) ,DU3 ( 3 ) , DFNB ,DFND, DFNF 
COMMON /DFN2/ AE , BE , DFNR , ZSTAR , AGL , BETA1 , BETA 2 , BETA3 , BETA4 
COMMON /DFN3/ AD, BD,BETA( 37 ) , DELRHO( 37 ) ,RADOC ( 37 ) 

COMMON /DFN4/ I HOLE, PL ( 6,6,4) , PA( 2,6,4) , EDG( 6,4) 

COMMON /DFN 5/ CE,CO,ABOT,ATOP 

COMMON /DFN 6/ RTADI , RTADO, RCI , RCO, RDI , RDO, XCI , XCO, XTADI , XT ADO 
COMMON /DFN7/ RW021 ,XW021 , RADII , RADIO 

C 

DIMENSION CADDO) ,VEC(3) ,VDUM(3) ,UE(3) 


GENERAL CONSTANTS 


I HOLE ■ 0 
— BOWL ELLIPSE 


AE - 4.5 
BE - 2.85776 


ORIGIN OF BOWL ELLIPSE 


DFND « 
DFNR - 

— BOWL ELLIPSE FOCI 


5.5 

6.02472 


FE ■ 

FIRST BOWL ELLIPSE 
XFE ■ 


SQRT(AE*AE - BE*BE ) 
FOCI AXIAL DISTANCE 
DFND - FE 


— ANGLES OF POINTS ON HOLE 

BETA1 - 130./RADDEG 
BETA2 - 50./RADDEG 
BETA3 - 298./RADDEG 
BETA4 - 220./RADDEG 

ORIGIN OF DUCT AXIS 

DFNB - 5.0 
DFNF ■ 5.2 


ANGLE OF DUCT 


AGL - 10.0/RADDEG 


UNIT VECTORS OF DUCT COORDINATE SYSTEM 


C 


DUl(l) - SIN(AGL) 
DUI ( 2 ) - 0 . 

DUI (3) » COS ( AGL ) 

DU2( 1 ) » COS (AGL) 
DU2 ( 2 ) » 0 . 

DU2 ( 3 ) » -SIN (AGL) 
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DU3(1) • 0. 

DU3 ( 2 ) » 1 . 

DU3{3) - 0. 

—DUCT ELLIPSE 

AD * 2.98 
BD » 3.55 

— AXIAL LENGTH TO DUCT EXIT PLANE . 

CE - 16.2 

— AXIAL LENGTH TO POINT ON EXIT PLANE NORMAL IN X 
CQ - 25.206 

— UNKNOWN 

ATOP - 7. 

ABOT - 3. 

— POINT 1: BOWL ENTRANCE - INNER 

RADII - 6.6 

— POINT 2: END OF BOWL - INNER 
XLI ■ 9.810 
RADOI - 6.029 

— POINT 3: END OF BOWL - OUTER 
XLO « 9.830 

CXL « ((XLO - DFND)/AE)**2 
RADOO - DFNR + BE*SQRT{1.0 - CXL) 

POINT 4: BOWL ENTRANCE - OUTER 




RADIO 

m 

7.5 

EDGE 

1: 

SEGMAX 1- 

mmmm 




XMAXll 

M 

2.450 

EDGE 

Is 

SEGMAX 3 





XMAX31 

m 

7.912 



RMAX31 

m 

6.380 

EDGE 

3: 

SEGMAX 1 





XMAX13 

m 

1.0 



RMAX13 

m 

7.815 



ZSTAR 

m 

RMAX13 
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— EDGE 3: SEGMAX 2 

CR23 - ( (RMAX13 - DFNR)/BE)**2 

XMAX23 - DFND - AE*SQRT(1.0 - CR23) 

EDGE 3: SEGMAX 3 

XMAX33 - 3.0 
XMAX33 « 2.35 

— EDGE 3: SEGMAX 4 

XMAX43 - DFND + AE*SQRT(1.0 - CR23) 

— EDGE 11; SEGMAX 3 

XMAX311 ■ 3.4 

CX311 - UXMAX3U - DFND)/AE)**2 
RMAX311 » DFNR + BE*SQRT(1.0 - CX311) 

EDGE 11: SEGMAX 4 

XMAX411 ■ 8.5 

CX411 » ((XMAX411 - DFND)/AE)**2 

RMAX411 - DFNR + BE*SQRT(1.0 - CX411) 

ANGLE OF INTERSECTION BETWEEN SECTION 1 AND SECTION 2 

THETAD - 75.0 

TH - THETAD/RADDEG 
CTH « COS(TH) 

STH - SIN(TH) 

— ANGLE OF INTERSECTION FOR SEGMENT 1 & 2 AND SEGMENT 2 ( 3 

THETAD 2 - 18.0 

C 

THE2 - THETAD2/RADDEG 
CTHE2 » COS(THE2) 

STHE2 ■ SIN(THE2) 

C 

THETAD3 - 82.0 

C 

THE3 ■ THETAD3/RADDEG 
CTHE3 ■ COS(THE3) 

STHE3 » SIN(THE3) 

C 

IF( IZONE.EQ. 1. AND. IZINDEX.GT. 1 ) THEN 
GO TO 40 
ELSE 

GO TO 50 
END IF 
C 
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40 GO TO (50,50,200,300,300,50) IZINDEX 

C 

50 IF(IZONE.GT.l) GO TO 100 
C 

c ******************************************* 
C ZONE 1 - (SIDE OF BOWL WITHOUT HOLE) 
c ******************************************* 

C EDGE COEFFICIENTS 

C 

C EDGE 3: ELLIPSE 

C 

COEFE( 1,3,3) ■ XFE 
COEFE( 3,3,3) ■ -DFNR 
COEFE(4,3,3) ■ 1.0 

— EDGE 5s CIRCULAR ARC 

COEFE (1,1,5) - 0.0 

COEFE( 4,1,5) » 1.0 

—EDGE 6: CIRCULAR ARC 

COEFE ( 1,1,6) - XLI 
COEFE (4, 1,6) ■ 1.0 

— EDGE 7: CIRCULAR ARC 


COEFE ( 1,1,7) ■ XLO 
COEFE (4, 1,7) ■ 1.0 

— EDGE 8: CIRCULAR ARC 

COEFE ( 1,1,8) - 0.0 

COEFE( 4,1,8) - 1.0 

— EDGE 11: ELLIPSE 

DO 10 K - 3,5 

COEFE (1,K,11) » XFE 
COEFE ( 2, K, 11) ■ DFNR*STH 
COEFE ( 3, K, 11) ■ -DFNR*CTH 
10 COEFE ( 4, K, 11) » 1.0 

— SURFACE COEFFICIENTS — 


COEFS(4,2) - 1.0 
COEFS(7,2) ■ 3.0 

COEFS(4,4) - 1.0 
COEFS( 7 ,4 ) ■ 3.0 

COEFS( 4 ,6 ) « 1.0 
COEFS( 7,6) « 3.0 


CORNER POINT COORDINATES 


POINT ( 1,1) 
POINT (2,1) 
POINT ( 3 , 1 ) 


0.0 

0.0 

-RADII 


l 
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C 


C 


C 


c 


G 


C 


c 


c 

c 

c 

c 


c 


c 



c 


c 


POINT( 1,2) 

m 

XLI 

POINT ( 2,2) 

m 

0.0 

POINT ( 3,2) 

m 

-RADOI 

POINT( 1,3) 

m 

XLO 

POINT( 2,3) 

* 

0.0 

POINT ( 3 , 3 ) 

m 

-RADOO 

POINT ( 1,4) 

m 

0.0 

POINT ( 2,4 ) 

m 

0.0 

POINT( 3,4) 

M 

-RADIO 

POINT ( 1,5) 

M 

0.0 

POINT ( 2,5) 

m 

RADII*STH 

POINT ( 3 , 5 ) 

s 

-RADII*CTH 

POINT ( 1,6) 

m 

XLI 

POINT ( 2*6) 

m 

RADOI*STH 

POINT( 3,6) 

m 

-RADOI *CTH 

POINT (1,7) 

m 

XLO 

POINT ( 2,7) 

m 

RADOO* STH 

POINT ( 3,7) 

m 

-RADOO*CTH 

POINT ( 1,8) 

m 

0.0 

POINT(2,8) 

m 

RADIO* STH 

POINT( 3*8) 

m 

-RADIO*CTH 


--EDGE SEGMENT COORDINATES 


SEGMAX( 1,1, 1) 

X 

XMAXll 

S EGMAX (3,1, 1) 

* 

-RADII 

SEGMAXt 1 » 2 , 1) 

m 

XMAX31 

SEGMAX( 3,2, 1) 

M 

-RMAX31 

S EGMAX ( 1,1, 3) 

m 

XMAX13 

SEGMAXt 3.1, 3) 

m 

-RMAX13 

SEGMAXt 1 » 2, 3) 

m 

XMAX23 

SEGMAXt 3,2, 3) 

m 

-RMAX13 

SEGMAXt l tit 9) 

m 

XMAXll 

SEGMAXt 2,1, 9) 

m 

RADII*STH 

SEGMAXt 3,1, 9) 

m 

-RADII*CTH 

SEGMAXt 1,2, 9) 

m 

XMAX31 

SEGMAXt 2,2, 9) 

m 

RMAX31*STH 

SEGMAXt 3, 2, 9) 

• 

-RMAX31*CTH 

SEGMAXt 1,1, 11) 

m 

XMAX13 

SEGMAXt 2, 1,11) 

• 

RMAX13*STH 

SEGMAXt 3, 1,11) 


-RMAX13*CTH 

SEGMAXt 1,2,11) 

m 

XMAX23 

SEGMAXt 2, 2, 11) 

m 

RMAX13*STH 

SEGMAXt 3, 2, 11) 

m 

-RMAX13*CTH 

SEGMAXt 1,3, 11) 

• 

XMAX311 
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SEGMAX( 2,3,11) 
SEGMAX( 3,3,11) 

SEGMAX( 1,4,11) 
SEGMAX( 2,4,11) 
SEGMAX( 3,4,11) 


RMAX311*STH 

-RMAX311*CTH 

XNAX411 

RMAX411*STH 

-RMAX411*CTH 


EDGE NODE DISTRIBUTION- 


ETAMAX ( 1 , 1) 
ETAMAX ( 2 , 1) 


ETAMAX(1, 
ETAMAX ( 2 , 

ETAMAX ( 1 , 
ETAMAX ( 2 , 


3) 

3) 

9) 

9) 


ETAMAX ( 1,11) 
ETAMAX (2,11) 
ETAMAX ( 3,11) 
ETAMAX (4,11) 


12.0 

44.0 

5.0 

10.0 

12.0 

44.0 

5.0 

10.0 

19.0 

48.0 


RETURN 


ZONE 2 - (SIDE OF BOWL WITH HOLE) 
100 IF ( IZONE .GT . 2 )GO TO 200 
—EDGE COEFFICIENTS 


EDGE 3: ELLIPSE 


DO 110 K ■ 3,5 


COEFE( 1 ,K, 3) 
COEFE ( 2 , K , 3) 
COEFE ( 3 , K , 3) 
110 COEFE ( 4 , K , 3) 

» XFE 

- DFNR*STH 

- -DFNR*CTH 

- 1.0 

EDGE 5: CIRCULAR 

ARC 

COEFE (1,1, 5) 
COEFE (4,1, 5) 

- 0.0 
■ 1.0 

COEFE (1,2, 5) 
COEFE ( 4,2, 5) 

- 0.0 
■ 1.0 

COEFE (1,3, 5) 
COEFE (4,3, 5) 

- 0.0 
- 1.0 

EDGE 6: CIRCULAR 

ARC 

COEFE( 1,1, 6) 
COEFE ( 4,1, 6) 

- XLI 
■ 1.0 

COEFE (1,2, 6) 
COEFE (4,2, 6) 

■ XLI 

■ 1.0 
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COEFE( 1,3, 6) ■ XLI 
COEFE (4,3, 6) - 1.0 

— EDGE 7: CIRCULAR ARC 

C0EFE(1,1, 7) « XLO 

COEFE (4,1, 7) - 1.0 

COEFE ( 1*2, 7) « XLO 

COEFE ( 4,2, 7) « 1.0 

COEFE ( 1,3, 7) * XLO 
COEFE ( 4 ,3# 7) - 1.0 

— EDGE 8s CIRCULAR ARC 

COEFE( 1,1, 8) ■ 0.0 
COEFE ( 4,1, 8) ■ 1.0 

COEFE(4,2, 8) - 1.0 
COEFE (4, 3, 8) - 1.0 

— EDGE 11; ELLIPSE 

COEFE (1,3,11) » 0.0 
COEFE( 2,3>11) « 0.0 
COEFE( 3 ,3*11) * 0.0 
COEFE (4,3,11) ■ 0.0 

COEFE ( 1,4*11) - 0.0 
COEFE (2,4,11) - 0.0 
COEFE ( 3 >4,11) ■ 0.0 
COEFE(4,4, 11) » 0.0 

COEFE (1,5,11) * XFE 
COEFE(2,5,ll) » 0.0 
COEFE(3,5,ll) - DFNR 
COEFE ( 4,5,11) ■ 1.0 

SURFACE COEFFICIENTS 


COEFS (4,2) « 1.0 
COEFS( 7 , 2 ) - 3.0 

COEFS (4,6) » 1.0 
COEFS (7,6) « 3.0 

— CORNER POINT COORDINATES' 



C 


POINT ( 1.1) 
POINT ( 2,1) 
POINT ( 3 , 1 ) 
POINT(4,l) 
POINT (5,1) 

POINT ( 1,2) 
POINT ( 2, 2) 
POINT (3, 2) 
POINT ( 4 , 2 ) 
POINT(5,2) 


0.0 

RADII*STH 

-RADII*CTH 

0.0 

0.0 


RADOI*STH 

-RADOI*CTH 

THETAF 

90.0 


C-70 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 



non o non noon 


LMSC-HEC TR F225968 


POINT (1,3) 


XLO 

POINT ( 2 , 3 ) 


RADOO*STH 

POINT( 3,3) 


-RADOO*CTH 

POINT( 4,3) 


THETAF 

POINT ( 5,3) 


90.0 

POINT( 1,4) 


0.0 

POINT ( 2, 4 ) 


RADIO* STH 

POINT( 3,4) 


-RADIO*CTH 

POINT( 4,4) 


12.831 

POINT( 5,4) 


- 7.491 

POINT (1,5) 


0.0 

POINT ( 2,5) 


0.0 

POINT ( 3,5) 


RADII 

POINT ( 4,5) 


0.0 

POINT (5,5) 


0.0 

POINT( 1,6) 


XL I 

POINT( 2,6 ) 


0.0 

POINT ( 3,6) 


RADOI 

POINT ( 4,6) 


0.0 

POINT ( 5,6) 


0.0 

POINT ( 1,7) 


XLO 

POINT ( 2,7) 


0.0 

POINT ( 3,7) 


RADOO 

POINT { 4,7) 


0.0 

POINT ( 5,7) 


1 

tn 

o 

• 

o 

POINT( 1,8) 


0.0 

POINT { 2,8) 


0.0 

POINT ( 3,8) 


RADIO 

POINT ( 4,8) 


0.0 

POINT( 5,8 ) 


14.7 


— UPPER LEFT CORNER ON HOLE — — 

CALL HOLE (BETA1, EDG (1,1) ,VDUM) 

RADI - SQRT(EDG( 2, 1 )**2 + EDG(3,1)**2) 
ANGl - ACOS(EDG( 2,1) /RADI) 

• — UPPER RIGHT CORNER ON HOLE 

CALL HOLE (BETA2,EDG( 1,2) ,VDUM) 

RAD2 ■ SORT ( EDG ( 2 , 2 ) * * 2 + EDG(3,2)**2) 
ANG2 - ACOS ( EDG (2,2) /RAD2 ) 


LOWER RIGHT CORNER ON HOLE 

CALL RHOS(BETA( 31) ,DELRHO( 31) ) 

CALL HOLE ( BETA3 , EDG ( 1,3) ,VDUM) 

RAD3 » SORT (EDG (2, 3)** 2 + EDG(3,3)**2) 
ANG3 « ACOS ( EDG (2,3) /RAD3 ) 
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— LOWER LEFT CORNER ON HOLE 

CALL RHOS ( BETA (23), DELRHO (23)) 

CALL HOLE (BETA4,EDG( 1,4) ,VDUM) 

RAD4 * SORT ( EDG ( 2 , 4 ) * * 2 + EDG(3,4)**2) 
ANG4 - ACOS ( EDG ( 2 , 4 ) /RAD4 ) 

— EDGE SEGMENT COORDINATES 


SEGMAX( 1,1, 

1) 

m 

XMAXll 

SEGMAX( 2,1, 

1) 

m 

RADII*STH 

SEGMAX (3,1, 

1) 

* 

-RADII*CTH 

SEGMAX( 1,2, 

1) 

m 

XMAX31 

SEGMAX ( 2,2, 

1) 

a 

RMAX31*STH 

SEGMAX ( 3,2, 

1) 

m 

-RMAX31*CTH 

SEGMAX( 1,1, 

3) 

m 

XMAX13 

SEGMAX (2,1, 

3) 

M 

RMAX13*STH 

SEGMAX( 3,1, 

3) 

m 

-RMAX13*CTH 

SEGMAX ( 1,2, 

3) 

m 

XMAX23 

SEGMAX ( 2,2, 

3) 

m 

RMAX13*STH 

SEGMAX (3, 2, 

3) 

u 

-RMAX13*CTH 

SEGMAX( 1,3, 

3) 

m 

XMAX311 

SEGMAX (2, 3, 

3) 

m 

RMAX311*STH 

SEGMAX (3, 3, 

3) 

m 

-RMAX3 1 1 *CTH 

SEGMAX( 1,4, 

3) 

m 

XMAX411 

SEGMAX ( 2,4, 

3) 

m 

RMAX411*STH 

SEGMAX (3,4, 

3) 

m 

-RMAX411*CTH 

SEGMAX (1,1, 

7) 

m 

XLO 

SEGMAX ( 2,1, 

7) 

m 

RADOO*CTHE2 

SEGMAX (3,1, 

7) 

m 

RADOO*STHE2 

SEGMAX ( 1,2, 

7) 

m 

XLO 

SEGMAX ( 2,2, 

7) 

m 

RADOO*COS ( ANG3 ) 

SEGMAX(3,2, 

7) 

m 

RADOO*SIN( ANG3 ) 

SEGMAX (1,1, 

8) 

m 

0.0 

SEGMAX ( 2,1, 

8) 

m 

RADIO*CTHE2 

SEGMAX( 3,1, 

8) 

m 

RADIO*STHE2 

SEGMAX (1,2, 

8) 

m 

0.0 

SEGMAX ( 2,2, 

8) 

u 

RADIO*CTHE3 

SEGMAX( 3,2, 

8) 

m 

RADIO*STHE3 

SEGMAX ( 1,1, 

9) 

m 

XMAXll 

SEGMAX( 2,1, 

9) 

m 

0.0 

SEGMAX (3,1, 

9) 

m 

RADII 

SEGMAX ( 1,2, 

9) 

9 

XMAX31 

SEGMAX ( 2,2, 

9) 

m 

0.0 

SEGMAX( 3 ,2, 

9) 

M 

RMAX31 
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SEGMAX( 1,1,11) 

m 

XMAX13 

SEGMAX ( 2,1,11) 

9 

o 

• 

o 

SEGMAX ( 3,1,11) 

m 

RMAX13 

SEGMAX ( 1,2,11) 

9 

XMAX23 

SEGMAX ( 2,2,11) 


0.0 

SEGMAX ( 3,2,11) 

9 

RMAX13 

SEGMAX( 1,3,11) 

9 

XMAX33 

SEGMAX ( 2,3,11) 

9 

0.0 

SEGMAX ( 3,3,11) 

9 

RMAX13 

SEGMAX (1,4, 11) 

9 

XMAX43 

SEGMAX ( 2,4,11) 

9 

0.0 

SEGMAX ( 3,4,11) 

9 

RMAX13 


DESCRIPTION OF COORDINATES ON SURFACE WITH HOLE (SIDE 4) 


DO 130 J » 1,3 
C 

PL(J,1,1) 
PL( J, 2, 1 ) 
PL(J,3,1) 
PL( J, 4 , 1 ) 
PL(J,5,1) 
PL( J , 6 , 1 ) 
C 

PL( J , 1 , 2 ) 
PL(J,4,2) 
PL( J ,5,2) 
PL(J,6,2) 
C 

PL( J ,1,3) 
PL( J , 4 , 3 ) 
PL( J ,5,3) 
PL( J,6,3) 
C 

PL( J ,1,4) 
PL( J, 2,4 ) 

PL( J ,3,4) 
PL( J, 4,4) 
PL(J ,5,4) 
130 PL(J,6,4) 

C 

PL( 1,2,2) 
PL( 2, 2, 2 ) 
PL(3,2,2) 

C 

PL( 1,3,2) 
PL (2, 3, 2) 
PL( 3 ,3,2) 

C 

PL( 1,2,3) 
PL( 2,2,3) 
PL(3,2,3) 
C 

PL(1,3,3) 
PL( 2,3,3) 
PL(3,3,3) 


POINT(J ,4) 
SEGMAX( J ,1,3) 
SEGMAX( J ,2,3) 
SEGMAX( J ,3,3) 
SEGMAX( J, 4,3) 
POINT (J , 3 ) 

SEGMAXt J ,1,8) 
EDG( J , 1 ) 

EDG( J , 2 ) 
SEGMAX( J ,1,7) 

SEGMAX( J ,2,8) 
EDG( J ,4) 

EDG( J ,3) 
SEGMAX( J ,2,7) 

POINT (J , 8 ) 
SEGMAX (J,l,ll) 
SEGMAX( J, 2, 11 ) 
SEGMAX (J ,3,11) 
SEGMAX( J , 4,11) 
POINT (J, 7) 

XMAX13 

RMAX13*CTHE2 

RMAX13*STHE2 

XMAX23 

RMAX13*CTHE2 

RMAX13*STHE2 

XMAX13 

RMAX13*CTHE3 

RMAX13*STHE3 

XMAX23 

RMAX13*CTHE3 

RMAX13*STHE3 
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C 

c 

C EDGE NODE DISTRIBUTION 


C 


c 


c 


c 

c 


c 

r 


ETAMAX ( 1 , 1) 

m 

12.0 

ETAHAX ( 2 , 1) 

m 

44.0 

ETAMAX ( 1 , 3) 

m 

5.0 

ETAMAX ( 2 , 3) 

m 

10.0 

ETAMAX ( 3 , 3) 

* 

19.0 

ETAMAX ( 4 , 3) 

m 

48.0 

ETAMAX ( 1 , 9) 

m 

12.0 

ETAMAX ( 2 , 9) 

m 

44.0 

ETAMAX (1,11) 

m 

5.0 

ETAMAX (2,11) 

m 

13.0 

ETAMAX( 3,11) 

m 

19.0 

ETAMAX (4,11) 

m 

48.0 

DO 140 K * 7,8 

ETAMAX ( 1 , K) 

m 

25.0 

ETAMAX ( 2 , K) 

• 

68.0 


RETURN 

c * ft******************************************' 

C 20NE 3 (DUCT) 

c «*******************************************< 

C 

200 CONTINUE 

GO TO (500,500,201,300,300,201) IZINDEX 

C 

201 CONTINUE 
C 

IF( I ZONE. GT. 3) GO TO 300 

C 

DO 205 1*1 , 8 
DO 205 J-1,6 
COEFS ( I , J ) « 0.0 
C 

DO 205 K-1,12 
205 COEFE(I,J,K) « 0.0 

C 

C CORNER POINT COORDINATES ON HOLE 


C 

C 


C 

C- 

c 

c 

c 


DO 210 J*l,3 

POINT (J , 1 ) » PL( J,4, 2) 

POINT ( J , 2 ) - PL( J, 5, 2) 

POINT( J , 6) - PL( J, 5,3) 

210 POINT( J, 5) - PL(J,4,3) 

—CORNER POINT COORDINATES ON EXIT PLANE 

CALL DEXIT(BETA1 , POINT ( 1,4) ,VDUM) 

CALL DEXIT(BETA2,POINT(l,3) ,VDUM) 
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CALL DEXIT(BETA3,P0INT( 1,7) ,VDUM) 

C 

CALL DEXIT(BETA4, POINT (1,8) ,VDUM) 
EDGE COEFFICIENTS ALONG HOLE (MAP « 7) 


COEFE( 1,1, 

1) 

m 

BETAl 

COEFE (2,1, 

1) 

m 

BETA2 

COEFE (1,1, 

5) 

m 

BETAl 

COEFE (2,1, 

5) 

m 

BETA4 

COEFE (1,1, 

6) 

m 

BETA2 

COEFE (2,1, 

<5) 

m 

BETA3 - 2.0*PI 

COEFE (1,1, 

9) 

m 

BETA4 

COEFE (2,1, 

9) 

m 

BETA3 


EDGE COEFFICIENTS 

COEFE (1,1, 2) 
COEFE (1,1, 4) 
COEFE( 1,1,10) 
COEFE( 1,1,12) 
EDGE COEFFICIENTS 


COEFE (1,1, 

3) 

m 

BETAl 

COEFE (2,1, 

3) 

m 

BETA2 

COEFE (1,1, 

7) 

m 

BETA2 

COEFE ( 2,1, 

7) 

m 

BETA3 

COEFE ( 1,1, 

8) 

m 

BETAl 

COEFE (2,1, 

8) 

m 

BETA4 

COEFE (1,1, 

11) 

m 

BETA4 

COEFE ( 2,1, 

11) 

m 

BETA3 


ALONG DUCT (MAP « 9) 

BETA2 
BETA1 
BETA 3 
BETA4 

ALONG EXIT ELLIPSE (MAP 


8 )- 


SEGMENT MAXIMUMS ALONG DUCT (CIRCULAR ARCS)- 

ISEG - 1 
RATIO ■ -1.0 

DO 250 J-1,4 

I EDGE • 4 
IF(J.EQ.2) IEDGE - 2 
IF( J.EQ.3) IEDGE - 10 
IF( J.EQ.4 ) IEDGE * 12 

EDG( 1 , J ) - COEFE (1, ISEG, IEDGE) 
MAPSEG( ISEG, IEDGE) • 9 


C 

C 


CALL EMAP ( IEDGE , ISEG , RATIO , EDG ( 1 , J ) , VDUM ) 

EDG( I ,J) 


DO 240 1-1,3 
240 SEGMAX( I, ISEG, IEDGE) 
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C 

MAXETA 
ETAMAX ( ISEG, IEDGE) 

C 

SEGMAX ( 4 , ISEG , IEDGE ) 
250 SEGMAX( 5, ISEG, IEDGE) 


0 . 65*NMBRNDS ( 2 ) 
MAXETA 

TXT 

TXZ 


— SURFACE COEFFICIENTS 


COEFS (1,1) 
COEFS (2,1) 
COEFS (3,1) 

COEFS (1,3) 
COEFS (2,3) 
COEFS (3, 3) 

COEFS (1,5) 
COEFS (2, 5) 
COEFS (3, 5) 

COEFS (1,6) 
COEFS ( 2,6) 
COEFS (3, 6) 

RETURN 


BETA1 

BETA2 

ETAMAX (1,2) 

BETA4 

BETA3 

ETAMAX (1,12) 

BETA1 

BETA4 

ETAMAX (1,4) 
BETA2 

BETA3 - 2.*PI 
ETAMAX ( 1,10) 


ZONE 465 (TURN AROUND DUCT) 


300 CONTINUE 

GO TO (500,500,500,301,301,301) IZINDEX 

301 CONTINUE 

IF( IZINDEX.GE . 4 .AND. IZONE. EQ. 2)GO TO 350 

I F ( I ZONE . GT • 4 ) GO TO 350 

— INITIALIZE INPUT FOR ZONE 3 TO ZERO- 

DO 325 N - 1,8 

DO 310 M « 1,5 
310 POINT (M,N) - 0.0 

DO 325 M - 1,6 
325 COEFS (N,M) - 0.0 

DO 340 N - 1,12 
DO 340 M ■ 1,5 

DO 330 J - 1,5 
330 SEGMAX(J,M,N) - 0.0 

DO 335 J - 1,8 
335 COEFE(J,M,N) ■ 0.0 

340 ETAMAX ( M , N ) « 0.0 
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C INITIALIZE CONSTANTS 

C 

C POINT Is ENTRANCE 

C 

350 XTAOI « -2.76 

C 

RDO - 5.6 

POINT 4: ENTRANCE 

RDI « 4.56 

INNER CIRCULAR ARC 

RTADI « 0.4570 
XCI - -4.331 
RCI - RDO + RTADI 
RWI21 » RDO + 2.0*RTADI 

OUTER CIRCULAR ARC 

RTADO - 1.36 

XCO - -4.35 
RCO - RDI + RTADO 


— TAD EXIT 


XTADO - 0.0 

— TANGENCY POINT ON OUTER SURFACE 
CALL TADWALL 


ANGLE 1 


THETAX » 15.0 

THX - THETAX/RADDEG 
CTH - COS (THX) 

STH - SIN(THX) 


— ANGLE 2 


CTH1 • COS (THE 2) 
STH1 - SIN(THE2) 


— ANGLE 3 


C 

C 


CTH 4 
STH4 

THETAC 

ANG1C 


COS (THE 3) 
SIN(THE3) 

90.0 - THETAX 

90.0 - THETAD2 
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C*’ 

C 

C*’ 

c 

c 

c 

c- 


ANG4C - 90.0 - THETAD3 
ZONE 4 (TURN AROUND DUCT) 

IF( IZINDEX.GE . 4 .AND. I ZONE . EQ. 2 )GO TO 400 
I F ( I ZONE .GT. 4 ) GO TO 400 
EDGE COEFFICIENTS 


C 

C EDGE 1: CIRCULAR ARC 

C 

COEFE (1,2, 1) - XCI 
COEFE( 3,2, 1) ■ -RCI 

C 

C EDGE 3: CIRCULAR ARC 

C 

COEFE (1,2, 3) ■ XCO 
COEFE (3, 2, 3) « -RCO 
C 

C EDGE 9s CIRCULAR ARC 


C 


c 

c 

c 


c 

c 


COEFE( 1,2, 9) ■ XCI 
COEFE ( 3 , 2 , 9) - RCI 

EDGE 11: CIRCULAR ARC 

COEFE(l,2,ll) * XCO 
COEFE ( 3*2*11) - RCO 

EDGE : CIRCULAR ARC 


C 

C 

c 

c- 


DO 360 J » 1,4 

COEFE (1,J, 5) - XT AD I 
COEFE ( 1 , J , 8) ■ XT AD I 
COEFE (1,J, 6) » XCI - 
360 COEFE (1,J, 7) « XCO - 

—SURFACE COEFFICIENTS — — 


RTADI 

RTADO 


C 

C 

c 


c 

c- 

c 



c 


COEFS(4,2) « 1.0 
COEFS (7,2) - 3.0 

COEFS (4,4) • 1.0 
COEFS (7, 4) » 3.0 

COEFS (4, 6) ■ 1.0 
COEFS (7,6) » 3.0 

—CORNER POINT COORDINATES 


POINT (1,1) » XT AD I 
POINT ( 2,1) ■ 0,0 

POINT( 3, 1 ) « -RDO 


POINT (1,2) ■ XCI - RTADI 
POINT ( 2,2) « 0,0 

POINT (3, 2) ■ -RCI 


C-78 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 


nnn o on onooo no no non 


LMSC-HEC TR F225968 


C 


SEGMAX( 2,2, 

7) - 

RCO*CTH 1 

SEGMAX( 3,2, 

7) - 

RC0*STH1 

SEGMAX( 2,3, 

7) - 

RCO*CTH4 

SEGMAX( 3,3, 

7) - 

RCO*STH4 

SEGMAX( 2,1, 

8) 

RDI*CTH 

SEGMAX( 3,1, 

8) 

-RDI*STH 

S EGMAX (2,2, 

8) 

RDI*CTH1 

SEGMAX( 3,2, 

8) 

RDI*STH1 

S EGMAX ( 2,3, 

8) 

RDI*CTH4 

SEGMAX( 3,3, 

8) 

RDI*STH4 

S EGMAX ( 1,1, 

9) 

XCI 

SEGMAX( 3,1, 

9) 

RDO 

S EGMAX (1,1, 

11) 

XCO 

SEGMAX( 3,1, 

11) 

RDI 


—EDGE NODE DISTRIBUTION — 

ETAMAXd, 1) - 13.0 

ETAMAXd, 3) ■ 13.0 

ETAMAXd, 9) « 13.0 

ETAMAXd, 11) - 13.0 

DO 380 J - 7,8 

ETAMAXd, J) ■ 36.0 
ETAMAX ( 2 , J) » 60.0 
380 ETAMAXd, J) « 103.0 

ETAMAXd, 5) - 36 
ETAMAXd, 6) - 36 

RETURN 


DATA FOR ZONE 5 (TURN AROUND DUCT) 


— INITIALIZE EDGE COEFFICIENTS 

400 COEFE(l , 2, 1) - 0.0 
COEFE( 3, 2, 1) ■ 0.0 

COEFE(l,2, 3) ■ 0.0 
COEFE(3,2, 3) ■ 0.0 

COEFE( 1,2, 9) ■ 0.0 
COEFE(2,2, 9) ■ 0.0 
C0EFE( 3, 2, 9) ■ 0.0 

C0EFE(1,2,11) ■ 0.0 
C0EFE( 2, 2, 11 ) « 0.0 
COEFE(3,2,ll) • 0.0 

— EDGE COEFFICIENTS 


PRECEDING PAGE BLANK NOT FILMED 
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EDGE Is CIRCULAR ARC 


COEFE (1,1, 1) • XCI 
COEFE (3,1, 1) - -RCI 


EDGE 3: CIRCULAR ARC 

COEFE (1,1, 3) ■ XCO 
COEFE (3,1, 3) - -RCO 

EDGE 9: CIRCULAR ARC 

COEFE (1,1, 9) - XCI 
COEFE (3,1, 9) - RCI 

— EDGE 11: CIRCULAR ARC 


COEFE (1,1,11) - XCO 
COEFE (3,1,11) - RCO 


— EDGE : CIRCULAR ARC 


DO 410 J - 1,4 


COEFE ( 1 , J , 5 ) 
COEFE ( 1 , J , 8 ) 
COEFE ( 1 , J, 6 ) 
410 COEFE ( 1 , J , 7 ) 


XCI - RTADI 
XCO - RTADO 
XTADO 
XT ADO 


— SURFACE COEFFICIENTS — 

COEFS (1,2) « 0.0 
COEFS (4,2) -1.0 
COEFS (7,2) - 3.0 

COEFS (1,4) - 0.0 
COEFS (4,4) ■ 1.0 
COEFS (7,4) - 3.0 

COEFS (1,5) - 0.0 
COEFS (4, 5) - 1.0 
COEFS (7,5) - 3.0 

COEFS (1,6) - 0.0 
COEFS (4,6) - 0.0 
COEFS (7,6) - 0.0 


CORNER NODE COORDINATES 


POINT( 1,1) 
POINT( 3 , 1 ) 
C 

POINT( 1,2) 
POINT ( 3,2) 
C 

POINT( 1,3) 
POINT ( 3,3) 
C 

POINT( 1,4) 
POINT ( 3,4) 


XCI - RTADI 
-RCI 

XTADO 

-RADII 

XTADO 

-RADIO 

XCO - RTADO 
-RCO 
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C 

c 


c 

c 


POINTU 

,5) 

• XCI - RTADI 

POINT (3 

,5) 

- RCI 

POINTU 

,6) 

* XTADO 

POINT (3 

,6) 

« RADII 

POINT (1 

,7) 

- XTADO 

POINT (3 

,7) 

« RADIO 

POINT (1 

,8) 

« XCO - RTADO 

POINT (3 

,8) 

« RCO 

)GE SEGMENT 

COORDINATES 

SEGMAX( 1,1, 

i) 

• XCI 

SEGMAX( 3,1, 

i) 

« -RWI21 

SEGMAX( 1,2, 

i) 

- -2.128 

SEGMAX( 3,2, 

i) 

« -6.514 

SEGMAX( 1,1, 

3) 

» XW021 

SEGMAX( 3,1, 

3) 

- -RW021 

SEGMAX( 1,2, 

3) 

« -2.2 

SEGMAX( 3,2, 

3) 

« -7.5 

DO 420 J ■ 

1,3 


SEGMAX( 1 , J , 

5) 

- XCI - RTADI 

SEGMAX ( 1 » J » 

8) 

* XCO - RTADO 

SEGMAX (2,1, 

5) 

■ RCI*CTH 

SEGMAX (3,1, 

5) 

« -RCI*STH 

SEGMAX ( 2,2, 

5) 

» RCI*CTH1 

SEGMAX ( 3,2, 

5) 

« RCI*STH1 

SEGMAX ( 2,3, 

5) 

- RCI*CTH4 

SEGMAX ( 3,3, 

5) 

« RCI*STH4 

DO 440 J « 

1,3 


SEGMAX ( 1 , J , 

6) 

■ 0.0 

SEGMAX ( 1 , J, 

7) 

■ 0.0 

SEGMAX(2,1, 

6) 

- RADII*CTH 

SEGMAX ( 3,1, 

6) 

• -RADII*STH 

S EGMAX (2,2, 

6) 

■ RADII*CTH1 

SEGMAX (3,2, 

6) 

• RADII*STH1 

SEGMAX (2, 3, 

6) 

■ RADII*CTH4 

SEGMAX (3, 3, 

6) 

- RADII*STH4 

SEGMAX ( 2,1, 

7) 

» RADIO*CTH 

SEGMAX ( 3 , 1 , 

7) 

» -RADIO*STH 

SEGMAX(2,2, 

7) 

- RADIO*CTHl 

SEGMAX ( 3,2, 

7) 

■ RADIO*STHl 
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SEGMAX( 2,3, 7) 

m 

RADIO*CTH4 

SEGMAXt 3,3, 7) 

m 

RADIO* STH4 

SEGMAXt 2,1, 8) 

m 

RCO*CTH 

SEGMAX( 3,1, 8) 

m 

-RCO*STH 

SEGMAX( 2,2, 8) 

M 

RCO*CTHl 

SEGMAX( 3,2, 8) 

m 

RCO*STHl 

SEGMAX( 2,3, 8) 

m 

RCO*CTH4 

S EGMAX (3,3, 8) 

m 

RCO*STH4 

SEGMAX( 1,1, 9) 

m 

XCI 

SEGMAX( 3,1, 9) 

m 

RWI21 

S EGMAX ( 1,2, 9) 

m 

•2.128 

S EGMAX ( 3,2, 9) 

• 

6.514 

S EGMAX ( 1,1,11) 

m 

XW021 

SEGMAX( 3,1,11) 

M 

RW021 

SEGMAXt 1,2,11) 

* 

•2.2 

SEGMAXt 3,2,11) 

m 

7.5 

—EDGE NODE DISTRIBUTION 

ETAMAX ( 1, 1) 

m 

13.0 

ETAMAX ( 1 , 3) 

m 

13.0 

ETAMAX ( 1 , 9) 

• 

13.0 

ETAMAX (1,11) 

m 

13.0 

ETAMAXt 2, 1) 

m 

31.0 

ETAMAX ( 2 , 9) 

m 

31.0 

ETAMAXt 2, 3) 

m 

31.0 

ETAMAXt 2, 11) 

m 

31.0 

DO 388 J • 7,8 

ETAMAXt 1, J) 

m 

36.0 

ETAMAXt 2, J) 

m 

60.0 

388 ETAMAXt 3, J) 

* 

103.0 


ETANAX (1,5) » 36 
ETAMAX (1,6) ■ 36 


500 CONTINUE 

RETURN 

END 


*****HGM 2*************************** 
********T**************************** 

SUBROUTINE HOLE (B, POINT, TANGENT) 


COMPUTE COORDINATES AND DERIVATIVE FOR A POINT ON THE HOLE. 
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B ANGULAR LOCATION OF A POINT ON THE HOLE 

POINT COORDINATES AT THE POINT 
TANGENT DERIVATIVE AT THE POINT 


COMMON /DFN1/ DU1 ( 3 ) ,DU2 ( 3 ) ,DU3 ( 3 ) ,DFNB,DFND,DFNF 
DIMENSION P0INT(6) »TANGENT(3) ,BC(3),US(3) 

—CALCULATE HOLE RADIUS 

CALL DELRAD ( B , RHO , DRHO , RHOD , RDOC ) 

— CALCULATE DUCT AXIAL DISTANCE 
CALL CAXIS ( B , RHO , C ) 

— ANGLULAR LOCATION OF A POINT ON THE HOLE 

SB - SIN(B) 

CB - COS(B) 

— UNIT VECTOR PERPENDICULAR FROM DUCT AXIS TO A POINT ON THE HOLE 

BC( 1 ) - CB*DU2 ( 1 ) + SB*DU3( 1 ) 

BC( 2 ) - CB*DU2 ( 2 ) + SB*DU3(2) 

BC ( 3 ) - CB*DU2 ( 3 ) + SB*DU3{3) 

VECTOR FROM BOWL CENTER TO A POINT ON 

XD « DFNF - DFND + C*DU1(1) + 

YD - DFNB + C*DU1(2) + 

ZD - C*DU1(3) + 

COORDINATES OF A POINT ON THE HOLE 

POINT ( 1 ) - XD + DFND 
POINT ( 2 ) - YD 
POINT ( 3 ) - ZD 

— UNIT VECTOR FROM BOWL CENTER TO A POINT ON THE HOLE 

PMAG - SQRT(XD*XD + YD* YD + ZD* ZD) 

US(1) - XD/PMAG 

US(2) - YD/PMAG 

US ( 3 ) - ZD/PMAG 

RETURN 
END 


THE HOLE 

RHO*BC( 1 ) 
RHO*BC ( 2 ) 
RHO*BC( 3 ) 


*****HGM 2***************************** 
********T ****************************** 

SUBROUTINE DUCT ( B , RHO , CAXIS , POINT ) 


COMPUTE DUCT POINT » X,Y,Z FOR GIVEN B,RHO,C 
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C 

COMMON /DFN1/ DU1( 3 ) , DU2 ( 3 ) , DU3 ( 3 ) ,DFNB ,DFND,DFNF 

C 

DIMENSION P0INT(6) ,BC(3) 


— POINT ANGLE 


CB « COS(B) 
SB » SIN(B) 


— RADIAL UNIT VECTOR 


BC ( 1 ) 
BC ( 2 ) 
BC ( 3 ) 


CB*DU2(1) 
CB*DU2 ( 2 ) 
CB*DU2 ( 3 ) 


+ SB*DU3 ( 1 ) 
+ SB*DU3 ( 2 ) 
+ SB*DU3( 3) 


— COORDINATES 


POINT(l) 
POINT ( 2 ) 
POINT( 3 ) 


DFNF + CAXIS*DUI(1) 
DFNB + CAXIS*DU1(2) 
CAXIS*DU1 ( 3 ) 


+ RHO*BC ( 1 ) 
+ RH0*BC(2) 
+ RHO*BC( 3) 


RETURN 

END 


HGM 2****************** 
*••▼*•*•*****«•***•**** 

SUBROUTINE RHOS(B,DRHO) 


COMPUTE THE DIFFERENCE IN RADIUS BETWEEN DUCT AND HOLE FOR THE 
LOWER CORNER POINTS ON THE HOLE. 

B ANGULAR LOCATION OF A POINT ON THE HOLE 

DRHO DIFFERENCE BETWEEN HOLE AND DUCT RADIUS 


COMMON /DFN1/ DU1( 3) ,DU2(3) ,DU3(3) , DFNB, DFND, DFNF 

COMMON /DFN2/ AE , BE , DFNR , ZSTAR , AGL . BETA1 , BETA2 , BETA3 , BETA4 

COMMON /DFN3/ AD, BD, BETA( 37 ) , DELRHO(37) ,RADOC(37) 

— TOLERANCE FOR NEWTON-RAPHSON ITERATION 

FEPS « 1.0E-07 


ANGLE OF DUCT AXIS 

SA - SIN (AGL) 

CA ■ COS (AGL) 

■ANGLULAR LOCATION OF A POINT ON THE HOLE 


SB * SIN(B) 

CB - COS ( B ) 

RADIUS OF ELLIPTIC DUCT 

FD = SQRT( (AD*SB)**2 + (BD*CB)**2) 
RHOD « AD*BD/FD 
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— INITIAL HOLE RADIUS 
RHO « RHOD ♦ DRHO 

— DISTANCE ALONG MAJOR AXIS OF BOWL ELLIPSE 

10 XMD ■ DFNF - DFND ♦ ZSTAR*SA/CA + RH0*CB/CA 
—ELLIPSE RATIO OF BOWL 
BOA * BE/AE 

—DISTANCE ALONG MINOR AXIS OF BOWL ELLIPSE 
FI » SQRT(AE**2 - XMD**2) 

— RADIAL DISTANCE TO HOLE 
F2 » DFNR + BOA* FI 
— VERTICAL DISTANCE TO HOLE 
F3 « DFNB + RHO* SB 

— HORIZONTAL DISTANCE TO HOLE - ZSTAR**2 
FR « F2**2 - F3**2 - ZSTAR**2 
DERIVATIVE 

DF1 « BOA*XMD*CB*F2/(Fl*CA) 

DF2 » F3*SB 
DFDR * -2. * ( DF1 + DF2) 

HOLE RADIUS 

RHO * RHO - FR/DFDR 

IF(ABS(FR) .GT.FEPS) GO TO 10 

DIFFERENCE IN RADIUS BETWEEN HOLE AND DUCT 

DRHO « RHO - RHOD 

RETURN 

END 

ft *******************#*************************** 1 

*****UTILITY******* ****************** ************ 

fr *******##*******##***************************** 1 

SUBROUTINE ANGLES ( VECTOR , THETA , PH I ) 

CONVERTS FROM VECTOR TO ANGLES 

THETA « THE ANGLE BETWEEN THE VECTOR 

AND ITS PROJECTION IN THE XZ PLANE 
PHI > THE ANGLE IN THE XZ PLANE 
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COMMON /INITA/ MAPTEN , INCHES 

COMMON /INITC/ PI, RADDEG 

C 

DIMENSION VECTOR ( 3 ) 

C 


C 


C 


C 

c 

c 


c 


c 

c 


c 


c 


DO 10 1-1,3 

10 IF(ABS(VECTOR(I) ).LT. 0.00001) VECTOR(I) - 0.0 

THETA - 0.0 
PHI - 0.0 

CALL VMAG( VECTOR, VECMAG) 

IF( VECMAG .EQ.0.0 ) RETURN 

20 VNORM2 - VECTOR ( 2 )/VECMAG 

IF( ABS ( VNORM2 ) . GT .1.0) VNORM2 - SIGN(ONE,VNORM2) 

THETA - AS I N ( VNORM 2 ) * RADDEG 
XZ PLANE • SQRT ( VECTOR ( 1 ) * * 2 + VECTOR(3)**2) 

IF ( XZ PLANE . LT. ABS (VECTOR ( 2))/1000.) RETURN 

VNORM 3 » VECTOR( 1 )/XZ PLANE 
IF (ABS ( VNORM3 ) .GT.1.0) VNORM3 - SIGN(ONE,VNORM3) 

PHI - ACOS(VNORM3)* RADDEG 
IF(VECTOR(3) . LT.0.0) PHI - -PHI 

RETURN 

END 


*****HGM 2 

********T* 


SUBROUTINE SPEDGE ( SRI , SRF , EPS , SAC , SEDGE ) 


CIRCULAR SUB-EDGE COORD., SIDE 4 OF BOWL 


DIMENS ION SRI ( 6 ) , SRF ( 6 ) , SEDGE ( 6 ) , SAC ( 3 ) , VDUM ( 3.) 
DIMENSION PCI (6) ,PC2(6) ,PVEC1(3) ,PVEC2(3) 
DIMENSION UNO) ,UP(3) ,UR(3) ,XN(3) ,XP(3),XR(3) 

C 

DATA PI /3. 141592654/ 

C 

DO 10 J-l , 6 
PC 1 ( J ) - SRI ( J) 

10 PC2( J) - SRF(J) 

C 

DO 20 J-l, 3 

PVECl(J) - PCI ( J ) - SACO) 

20 PVEC2(J) ■ PC2( J) - SACO) 

C 

CALL VMAG( PVEC1 , RM1 ) 

C 
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CALL VMAG(PVEC2,RM2) 

CALL CROSS ( PVEC 2 , PVEC 1 , UN , 7 1 ) 

CALL CR0SS(PVEC1,UN,UP,72) 

CALL CROSS (UN, UP, UR, 73) 

RCC ■ RM1*RM2 

CALL VDOT { PVEC i , PVEC 2 , RR ) 

THETA - PI 

IF(ABS(RR) . LE. RCC) THETA - ACOS(RR/RCC) 

RC ■ I.O/SIN(THETA) 

THETA1 « THETA 

ANG1 * (1.0 - EPS) ‘THETA 
ANG2 * EPS ‘THETA 

CANG1 - COS(ANGl) 

SANG1 » SIN(ANGl) 

CANG2 * COS( ANG2) 

SANG2 » SIN( ANG2 ) 

DO 40 J*1 , 3 

40 SEDGE(J) - SAC(J) + RC*(SANG1*PVEC1( J) + SANG2‘PVEC2( J) ) 

CALL VDOT (PCI (4) ,PC2(4) ,RR) 

IF(RR.GT. 0.9999) GO TO 70 

CALL CROSS(PC2(4) ,PC1 ( 4 ) ,XN ,74 ) 

CALL CROSS(PCl(4) ,XN,XP,75) 

CALL CROSS ( XN , XP , XR ,76) 

THET * ACOS(RR) 

ALPH » EPS*THET 

CANG • COS (ALPH) 

SANG ■ SIN (ALPH) 

CALL VADD ( CANG , XR , SANG , XP , VDUM , SEDGE ( 4 ) ) 

RETURN 
70 CONTINUE 
DO 80 J»l,3 

80 SEDGE(J+3) * PCllJ+3) 

RETURN 

END 
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C 

c*****«*t ******** ******* 

C*****hgm 2************* 

C ****....T*. ************ 

C 

SUBROUTINE TADWALL 


COMPUTE TANGENT POINT ON OUTER TAD WALL 


C 


C 


C 


C 

C 


C 


C 


COMMON /DFN6/ RTADI , RTADO , RCI , RCO , RDI , RDO , XCI , XCO , XTADI , XTADO 
COMMON /DPN7/ RW021 ,XW021 , RADII , RADIO 

TASQ » (XTADO - XCO)**2 + (RADIO - RCO)**2 - RTADO**2 
TAL » SORT ( TASQ ) 

ETA « TAL/RTADO 
ETASQ « ETA* *2 

A - 1.0 + ETASQ 
B - RADIO + ETASQ* RCO 
BSQ - B**2 

C * TASQ - RADIO**2 - ETASQ*RCO**2 
DSQ » 1.0 + (A*C/BSQ) 

D - SORT (DSQ) 

RW021 - B* ( 1 .0 + D)/A 

XW021 « XTADO - (RW021 - RCO)*ETA 

RETURN 

END 
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Appendix D 

HGM OUTPUT REWRITE LISTING 
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C 


SUBROUTINE REWRITE 


THIS PROGRAM CONVERTS HGM2DUCT CODE GEOMETRY FILES TO 
PL0T3D GEOMETRY FILES 


DIMENSION IDIM( 110) , JDIM( 110 ) ,KDIM( 110 ) 

DIMENSION XBUF( 200000) , YBUF ( 200000 ) ,ZBUF( 200000 ) 


MATE « 0 
— ASSIGN UNITS 

DATA 120, INPUT, ISCR,IP3D /1,2,3,4/ 

OPEN (UNIT* I 20, FILE* 'HGM2DUCT.DAT' jSTATUS-'OLD* ) 

OPEN ( UN IT* INPUT , FILE* * REWRITE .OUT 1 , STATUS- ' NEW’ ) 

OPEN ( UN IT* ISCR, STATUS* ' SCRATCH ’ ,FORM*' UNFORMATTED' ) 

— COUNTERS 

12 * 0 
NGRID ■ 0 

READ PARAMETERS FROM GEOMETRY 

I ERR * 1 

READ( 120,1000, END-200 ) NSTORE , IPLN, 

& INOD2 , JNOD2 , KNOD2 , MARCH2 

— -WRITE DATA RANGE 

11 - 12 + 1 

12 » 12 + NSTORE 

WRITE (INPUT, 1210) 

WRITE ( INPUT, 1220 ) NSTORE , IDYN, IPLN , INOD2 , JNOD2 , KNOD2 , 
& MARCH 2 , 11,12 


C 

C 

c 

c 

c 

c 

c 

c 


READ X,Y « Z FROM GEOMETRY 
I ERR * 2 

READ( 120,1010, ERR* 400 ) (XBUF(I) ,1*11,12) 
I ERR - 3 

READ( 120,1010, ERR-400 ) ( YBUF ( I ) , I-Il , 12 ) 
I ERR - 4 

READ( 120, 1010, ERR-400) (ZBUF(I) ,1-11,12) 
I ERR - 5 
I PLANE - 1 
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READ NEXT SET OF PARAMETERS FROM GEOMETRY 


100 READ( 120, 1000, END-200) NSTORE,IPLN,INOD,JNOD,KNOD, MATE, MARCH 
—CHECK FOR NEW ZONE 


C 


C 

C 

C 


c 


IF ( INOD.NE . INOD2 .OR. JNOD.NE. JNOD2 .OR. KNOD . NE . KNOD2 ) THEN 


I F ( MARCH 2 . EQ . 1 .AND. 
I F ( MARCH 2 . EQ . 2 .AND. 
I F ( MARCH 2 . EQ . 3 .AND. 


JNOD.EQ.JNOD2 .AND. KNOD. EQ. KNOD 2) 
INOD.EQ.INOD2 .AND. KNOD . EO « KNOD2 ) 
JNOD.EQ. JNOD2 .AND. INOD.EQ. INOD2) 


GO TO 110 
GO TO 110 
GO TO 110 


NGRID - NGRID + 1 
I F ( MARCH 2 . EO . 1 ) THEN 

I NOD 2 - IPLANE 
IDIM( NGRID) - KNOD2 
JDIM( NGRID) ■ JNOD2 
KDIM( NGRID) » INOD2 

END IF 


C 

I F ( MARCH 2 • EQ . 2 ) THEN 
C 

JNOD2 » IPLANE 
IDIM(NGRID) » INOD2 
JDIM( NGRID) - KNOD2 
KDIM( NGRID) - JNOD2 

C 

END IF 

C 

I F ( MARCH 2 . EO « 3 ) THEN 
C 

KNOD2 - IPLANE 
IDIM(NGRID) - JNOD2 
JDIM( NGRID) « INOD2 
KDIM( NGRID) ■ KNOD2 


C 

WRITE( INPUT, 1230) 

WRITE ( INPUT ,1220) NGRI D , IDIM ( NGRID ) , JDIM ( NGRID ) , KDIM ( NGRID ) 
C 

IERR - 6 
C 

WRITE (ISCR,ERR-411) 

1 (XBUF(I), 1-1,12), 

2 (YBUF(I), 1-1,12), 

3 (ZBUF(I) ,1-1,12) 

C 

WRITE( INPUT, 1250) 

C 

INOD2 - INOD 
JNOD2 - JNOD 
KNOD2 - KNOD 
C 

MARCH 2 - MARCH 
C 

WRITE ( INPUT, 1210) 
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C 


12 » 0 
IPLANE - 0 


END IF 


WRITE DATA RANGE 

110 CONTINUE 

11 - 12 + 1 

12 « 12 + NSTORE 

WRITE ( I NPUT , i 2 20 ) NSTORE , IDYN, IPLN , INOD , JNOD , KNOD , MARCH ,11,12 

READ X, Y & Z FROM GEOMETRY 

I ERR » 7 

READ( 120, 1010, ERR-400) (XBUF( I ) , 1*11 , 12 ) 

I ERR - 8 

READ( 120, 1010, ERR-400) ( YBUF( I ) , 1*11 , 12) 

I ERR « 1 

READ( 120, 1010, ERR-400) ( ZBUF( I ) , I-Il , 12 ) 

IPLANE - IPLANE + 1 
GO TO 100 

END OF GEOMETRY DATA 


200 
C 

C 

C 


C 

C 

C 


C 

C 

C 


CONTINUE 


NGRID - NGRID + 1 

IF( MARCH 2 . £Q» 1 ) THEN 

IN0D2 « IPLANE 
IDIM( NGRID) - KN0D2 
JDIM( NGRID) * JN0D2 
KDIM( NGRID) - IN0D2 

END IF 

I F ( MARCH 2 . EQ • 2 ) THEN 

JN0D2 • IPLANE 
IDIM(NGRID) ■ IN0D2 
JDIM(NGRID) ■ KN0D2 
KDIM( NGRID) > JN0D2 

END IF 

IF(MARCH2.EQ.3) THEN 

KN0D2 > IPLANE 
IDIM( NGRID) « JN0D2 
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C 


JDIM(NGRID) » IN0D2 
KDIM(NGRID) - RN0D2 

END IF 


WRITE ( INPUT, 1230) 

WRITE( INPUT, 1220) NGRID, IDIM(NGRID) , JDIM(NGRID) , KDIM(NGRID) 
WRITE (INPUT, 1240) 


— WRITE X , Y & 2 ON SCRATCH FILE 
IERR - 10 

WRITE ( ISCR,ERR-411 ) 

1 ( XBUF ( I ) , I-l , 12 ) , 

2 (YBUF(I) ,1-1,12) , 

3 (ZBUF(I), 1-1,12) 

WRITE (INPUT, 1250) 

CLOSE ( UN IT* I 20 ) 

WRITE( INPUT, 1260) 

REWIND( UNIT-ISCR) 

WRITE( INPUT, 1270) 


WRITE BINARY FILE 

OPEN (UN IT- IP 3D, FILE- ' PLT3D.BIN' , STATUS- ‘NEW* , FORM- 'UNFORMATTED 1 ) 
WRITE ( INPUT , 1300 ) 

IF(NGRID.GT.l) THEN 

IERR - 11 

WRITE( INPUT, 1310) NGRID 
WRITE( IP3D, ERR-420 ) NGRID 

END IF 


IERR - 11 
WRITE (INPUT, 1320) 

WRITE( INPUT, 1330) ( IDIM(N) , JDIM(N) , KDIM(N) ,N-1 , NGRID) 

WRITE (IP3D, ERR-420) ( IDIM(N) , JDIM(N) ,KDIM(N) ,N-1 , NGRID) 


— WRITE BINARY PLT3D FILE 


DO 300 N-l, NGRID 
C 

IERR - 12 
C 

12 « IDIM( N) *JDIM( N) *KDIM(N) 
C 
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READ(ISCR,ERR-410) 

1 

2 

3 

C 

IERR - 13 

C 

WRITE(IP3D,ERR-420) 

1 

2 

3 

C 

WRITE ( INPUT, 1340) N 

C 

300 CONTINUE 
—WRITE GRID NUMBER 


(XBUF(I) ,1-1,12), 
(YBUF(I) ,1-1,12) , 
(ZBUF(I) ,1-1,12) 


(XBUF(I), 1-1,12), 
(YBUF(I), 1-1,12), 
(ZBUF(I) ,1-1,12) 


WRITE( INPUT, 1350) NGRID 


GO TO 500 


ERRORS 


— ERROR READING GEOMETRY DATA FROM FILE20 
400 WRITE( INPUT, 1400) IERR 
GO TO 500 

— ERROR READING SCRATCH FILE 

410 WRITE( INPUT, 1410) IERR 

411 WRITE( INPUT, 1411) IERR 

GO TO 500 

— ERROR WRITING TO PLOT3D BINARY FILE 
420 WRITE( INPUT, 1420) IERR 


END OF PROGRAM 


500 CONTINUE 

CLOSE (UNIT-ISCR) 

CLOSE ( UNIT- IP 3D) 

CLOSE ( UNIT-INPUT ) 

FORMAT STATEMENTS 

1000 FORMAT (91 5) 

1010 FORMAT( 6E22. 14) 

C 

1200 FORMAT ( / 1 *** ERROR: RERUN GEOMN WITH MATE SET TO 0 ***') 
1210 FORMAT (/'NSTORE IDYN IPLN INOD JNOD KNOD MARCH II 12') 
1220 FORMAT ( IX, 715, 216) 

1230 FORMATt/' NGRID IDIM JDIM KDIM ') 

1240 FORMATt/' END OF FILE REACHED ON 120') 

1250 FORMATt/' GRID WRITTEN TO ISCR ') 
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1260 FORMATt/' 120 CLOSED') 
1270 FORMATt/' ISCR REWIND’) 


C 

1300 FORMATt/' IP3D OPENED') 

1310 FORMATt/' NGRID •',13) 

1320 FORMATt/' IDIM JDIM KDIM ') 

1330 FORMATt/ 315) 

1340 FORMATt/' GRID ',12,' WRITTEN TO PLT3D.BIN' ) 

1350 FORMATt/' THERE ARE ',12,' GRIDS WRITTEN TO PLT3D.BIN ) 

c 

1400 FORMATt/' ERROR READING GEOMETRY DATA FROM FILE20 ON UNIT 120', 12) 

1410 FORMATt/' ERROR READING SCRATCH FILE ON UNIT ISCR', 12) 

1411 FORMATt/' ERROR WRITING SCRATCH FILE TO UNIT ISCR' ,12) 

1420 FORMATt/' ERROR WRITING TO PLOT3D BINARY FILE ON UNIT IP3D',I2) 


STOP 

END 


C 
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